#include "normal_order.h"

int phase(int n) { return n % 2 ? -1 : 1; }

int Norm_Ord::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}
int Norm_Ord_Sec::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

void Norm_Ord::build_T2P_Tz(int a_beg, int a_end, int b_beg, int b_end) {
  // Emax = Emax_t;
  int j2a_max_a = 1;
  int j2a_max_b = 1;

  for (int a = a_beg; a <= a_end; a++) {
    int j2_t = SP.state_Tz[a].j_2;
    if (j2_t > j2a_max_a)
      j2a_max_a = j2_t;
  }
  for (int b = b_beg; b <= b_end; b++) {
    int j2_t = SP.state_Tz[b].j_2;
    if (j2_t > j2a_max_b)
      j2a_max_b = j2_t;
  }

  int J2_max = j2a_max_a + j2a_max_b;
  int sp_size = SP.state_Tz.size();
  int chan_size = (J2_max + 2) / 2;
  chan_size = 3 * chan_size; // pp , pn, nn Tz = 1, 0, -1
  state_T2P_Tz.resize(chan_size);

  cout << "Test size = " << sp_size << "\t ch_size = " << chan_size << endl;
  cout << a_beg << "\t" << a_end << "\t\t" << b_beg << "\t" << b_end << endl;

  State_DP state_t;

  int chan_num = -1;
  for (int Tz_2 = -2; Tz_2 <= 2; Tz_2 = Tz_2 + 2) {
    for (int J_2 = 0; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      for (int a = a_beg; a <= a_end; a++) {
        for (int b = b_beg; b <= b_end; b++) {
          if (Tz_2 != (SP.state_Tz[a].tz_2 + SP.state_Tz[b].tz_2))
            continue;
          // cout << "\t 2 : Tz_2 = " << Tz_2 << "\t J_2 = " << J_2
          //     << "\t a = " << a << "\t b = " << b << endl;
          // cout<<"Tri = "<<Tri(J_2, SP.state_Tz[a].j_2,
          // SP.state_Tz[b].j_2)<<endl;
          if (Tri(J_2, SP.state_Tz[a].j_2, SP.state_Tz[b].j_2) != 1)
            continue;
          // cout << "\t 3: Tz_2 = " << Tz_2 << "\t J_2 = " << J_2
          //     << "\t a = " << a << "\t b = " << b << endl;
          // Antisymmetry condition (1+(-1)^{Jab}\delta_{ab}) != 0
          int ph_Jab = J_2 / 2;
          if (a == b && (ph_Jab % 2 == 1))
            continue;
          state_t.a = a;
          state_t.b = b;
          state_t.J_2 = J_2;
          state_t.T_2z = Tz_2;
          state_t.E = SP.state_Tz[a].E + SP.state_Tz[b].E;
          state_T2P_Tz[chan_num].push_back(state_t);
        }
      }
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_T2P_Tz[i].begin(), state_T2P_Tz[i].end(),
                [](State_DP a, State_DP b) {
                  return (100000 * a.E + 100 * a.a + a.b) <
                         (100000 * b.E + 100 * b.a + b.b);
                });
    int index = 0;
    for (int j = 0; j < state_T2P_Tz[i].size(); j++) {
      state_T2P_Tz[i][j].index = index;
      index++;
    }
  }

  channel_T2P_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_T2P_Tz[i].size() == 0) {
      channel_T2P_Tz[i].num = 0;
      channel_T2P_Tz[i].J_2 = 0;
      channel_T2P_Tz[i].T_2z = -10;
      channel_T2P_Tz[i].E_min = 0;
      channel_T2P_Tz[i].E_max = 0;
    } else {
      channel_T2P_Tz[i].num = state_T2P_Tz[i].size();
      channel_T2P_Tz[i].J_2 = state_T2P_Tz[i][0].J_2;
      channel_T2P_Tz[i].T_2z = state_T2P_Tz[i][0].T_2z;
      int E_min = state_T2P_Tz[i][0].E;
      int E_max = state_T2P_Tz[i][0].E;
      for (int j = 0; j < state_T2P_Tz[i].size(); j++) {
        if (E_min > state_T2P_Tz[i][j].E)
          E_min = state_T2P_Tz[i][j].E;
        if (E_max < state_T2P_Tz[i][j].E)
          E_max = state_T2P_Tz[i][j].E;
      }
      channel_T2P_Tz[i].E_min = E_min;
      channel_T2P_Tz[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel_T2P_Tz[i].E_min;
    int Emax = channel_T2P_Tz[i].E_max;
    channel_T2P_Tz[i].indexE.resize(Emax + 1);
    channel_T2P_Tz[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state_T2P_Tz[i].size(); j++) {
      channel_T2P_Tz[i].indexE[state_T2P_Tz[i][j].E] = state_T2P_Tz[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state_T2P_Tz[i].size(); j++) {
        if (state_T2P_Tz[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel_T2P_Tz[i].indexE_beg[E] = indexE_beg;
    }
  }
}

int Norm_Ord::FindC_T2P_Tz(int J2ab, int T2ab_z) {
  int res = -1;
  int size_ch = channel_T2P_Tz.size();
  for (int ch = 0; ch < size_ch; ch++) {
    if (channel_T2P_Tz[ch].J_2 == J2ab && channel_T2P_Tz[ch].T_2z == T2ab_z) {
      res = ch;
      break;
    }
  }
  if (res < 0) {
    cout << "Wrong happened at FindC_T2P_Tz" << endl;
    cout << "J2ab = " << J2ab << "\t T2ab_z = " << T2ab_z << endl;
    exit(0);
  }
  return res;
}
int Norm_Ord::Find_T2P_Tz(int J2ab, int T2ab_z, int a, int b) {
  int ch_Tz = FindC_T2P_Tz(J2ab, T2ab_z);
  int size = state_T2P_Tz[ch_Tz].size();
  int index = -1;
  for (int i = 0; i < size; i++) {
    int a_t = state_T2P_Tz[ch_Tz][i].a;
    int b_t = state_T2P_Tz[ch_Tz][i].b;
    if (a == a_t && b == b_t) {
      index = i;
      break;
    }
  }
  if (index < 0) {
    cout << "Wrong happened at Find_T2P_Tz" << endl;
    cout << "J2ab = " << J2ab << "\t T2ab_z = " << T2ab_z
         << "\t ch_Tz = " << ch_Tz << "\t a = " << a << "\t b = " << b << endl;
    exit(0);
  }
  return index;
}

double Norm_Ord::cal1f3_Tz_C(State_SP &sp_t, int h_core) {
  double res = 0.0;

  int j2_p = sp_t.j_2;
  int t2z_p = sp_t.tz_2;
  int p = sp_t.index;

  for (int h1 = 0; h1 <= h_core; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    for (int h2 = 0; h2 <= h_core; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      int J2ab_min = std::abs(j2_h1 - j2_h2);
      int J2ab_max = j2_h1 + j2_h2;
      int T2z = t2z_p + t2z_h1 + t2z_h2;
      for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
        int J2_min = std::abs(J2ab - j2_p);
        int J2_max = J2ab + j2_p;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          int ch_Tz = TP.FindC_Tz(J2, T2z);

          int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, p);
          int y = x; // TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
          if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
            cout << "Wrong happened @ Norm_Ord::cal1f3_Tz" << endl;
            cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
            exit(0);
          }
          // double mat_t = this->vec_Tz[ch_Tz][x][y];
          double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
          double C1 = (J2 + 1.0) / (j2_p + 1.0);
          double mat = C1 * mat_t;

          double N_fac = 1;
          /*if (h1 == h2)
            N_fac = 0.5;/**/

          // cout << "\t J2 = " << J2 << "\t T2z = " << T2z
          //     << "\t J2ab = " << J2ab<<"\t j2_p = "<<j2_p<<"\t";
          //   << endl;
          // cout << "h1 = " << h1 << "\t h2 = " << h2 << "\t p = " << p
          // cout<<"C1 = "<<C1 << "\t mat_t = " << mat_t <<"\t mat = "<<
          // mat << endl; /**/
          res += mat * N_fac;
        }
      }
    }
  }
  res = 0.5 * res;
  return res;
}
double Norm_Ord::cal1f3_Tz_P(State_SP &sp_t, int h_P_beg, int h_P_end,
                             vector<double> rho_p_t, vector<double> rho_n_t) {
  double res = 0.0;

  int j2_p = sp_t.j_2;
  int t2z_p = sp_t.tz_2;
  int p = sp_t.index;
  double rho_p, rho_n;
  int flag = -1;
  if (t2z_p == -1) {
    // t2z = -1 n
    rho_n = rho_n_t[0];
    rho_p = rho_p_t[0];
    flag = 1;
  } else {
    // t2z = 1 p
    rho_n = rho_n_t[1];
    rho_p = rho_p_t[1];
    flag = 2;
  }
  /*cout<<"\t ====== rho_n = "<<rho_n<<"\t  rho_p = "<<rho_p<<"\t flag =
  "<<flag<<endl;
  cout<<rho_n_t[0]<<"\t"<<rho_n_t[1]<<"\t"<<rho_p_t[0]<<"\t"<<rho_p_t[1]<<endl;/**/
  int num_sum_p1, num_sum_n1;
  int num_sum_p2, num_sum_n2;

  for (int h1 = h_P_beg; h1 <= h_P_end; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    num_sum_n1 = 0;
    num_sum_p1 = 0;
    double rho_h1 = 0.0;
    if (t2z_h1 == -1) {
      rho_h1 = rho_n;
      num_sum_n1 = 1;
    } else {
      rho_h1 = rho_p;
      num_sum_p1 = 1;
    }
    for (int h2 = h_P_beg; h2 <= h_P_end; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      num_sum_p2 = 0;
      num_sum_n2 = 0;
      double rho_h2 = 0.0;
      if (t2z_h2 == -1) {
        rho_h2 = rho_n;
        num_sum_n2 = 1;
      } else {
        rho_h2 = rho_p;
        num_sum_p2 = 1;
      }
      int sum_n = num_sum_n1 + num_sum_n2;
      int sum_p = num_sum_p1 + num_sum_p2;

      if (sum_n > Num_valence_n)
        continue;
      if (sum_p > Num_valence_p)
        continue;

      int J2ab_min = std::abs(j2_h1 - j2_h2);
      int J2ab_max = j2_h1 + j2_h2;
      int T2z = t2z_p + t2z_h1 + t2z_h2;
      for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
        int J2_min = std::abs(J2ab - j2_p);
        int J2_max = J2ab + j2_p;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          int ch_Tz = TP.FindC_Tz(J2, T2z);

          int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, p);
          int y = x; // TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
          if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
            cout << "Wrong happened @ Norm_Ord::cal_val1f3_Tz 2" << endl;
            cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
            exit(0);
          }
          // double mat_t = this->vec_Tz[ch_Tz][x][y];
          double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
          double C1 = (J2 + 1.0) / (j2_p + 1.0);
          double mat = C1 * mat_t * rho_h1 * rho_h2;

          // cout << "\t J2 = " << J2 << "\t T2z = " << T2z
          //     << "\t J2ab = " << J2ab<<"\t j2_p = "<<j2_p<<"\t";
          //   << endl;
          /*cout << "h1 = " << h1 << "\t h2 = " << h2 << "\t p = " <<
          p<<endl; cout << "C1 = " << C1 << "\t mat_t = " << mat_t <<"\t
          rho1 =
          "<<rho_h1<<"\t rho2 = "<<rho_h2<<endl; cout<< "\t mat = " <<
          mat<< endl; /**/
          res += mat;
        }
      }
    }
  }
  res = 0.5 * res;
  return res;
}
double Norm_Ord::cal1f3_Tz_1C1P(State_SP &sp_t, int h_core, int h_P_beg,
                                int h_P_end, vector<double> rho_p_t,
                                vector<double> rho_n_t) {
  double res = 0.0;

  int j2_p = sp_t.j_2;
  int t2z_p = sp_t.tz_2;
  int p = sp_t.index;
  double rho_p, rho_n;
  int flag = -1;
  if (t2z_p == -1) {
    // t2z = -1 n
    rho_n = rho_n_t[0];
    rho_p = rho_p_t[0];
    flag = 1;
  } else {
    // t2z = 1 p
    rho_n = rho_n_t[1];
    rho_p = rho_p_t[1];
    flag = 2;
  }
  /*cout<<"\t ====== rho_n = "<<rho_n<<"\t  rho_p = "<<rho_p<<"\t flag =
  "<<flag<<endl;
  cout<<rho_n_t[0]<<"\t"<<rho_n_t[1]<<"\t"<<rho_p_t[0]<<"\t"<<rho_p_t[1]<<endl;/**/
  double res1 = 0.0;
  for (int h1 = h_P_beg; h1 <= h_P_end; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    double rho_h1 = 0.0;
    if (t2z_h1 == -1) {
      rho_h1 = rho_n;
    } else {
      rho_h1 = rho_p;
    }
    for (int h2 = 0; h2 <= h_core; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;

      int J2ab_min = std::abs(j2_h1 - j2_h2);
      int J2ab_max = j2_h1 + j2_h2;
      int T2z = t2z_p + t2z_h1 + t2z_h2;
      for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
        int J2_min = std::abs(J2ab - j2_p);
        int J2_max = J2ab + j2_p;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          int ch_Tz = TP.FindC_Tz(J2, T2z);

          int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, p);
          int y = x; // TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
          if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
            cout << "Wrong happened @ Norm_Ord::cal_val1f3_Tz 2" << endl;
            cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
            exit(0);
          }
          // double mat_t = this->vec_Tz[ch_Tz][x][y];
          double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
          double C1 = (J2 + 1.0) / (j2_p + 1.0);
          double mat = C1 * mat_t * rho_h1;

          // cout << "\t J2 = " << J2 << "\t T2z = " << T2z
          //     << "\t J2ab = " << J2ab<<"\t j2_p = "<<j2_p<<"\t";
          //   << endl;
          /*cout << "h1 = " << h1 << "\t h2 = " << h2 << "\t p = " <<
          p<<endl; cout << "C1 = " << C1 << "\t mat_t = " << mat_t <<"\t
          rho1 =
          "<<rho_h1<<"\t rho2 = "<<rho_h2<<endl; cout<< "\t mat = " <<
          mat<< endl; /**/
          res1 += mat;
        }
      }
    }
  }
  double res2 = 0.0;
  for (int h1 = 0; h1 <= h_core; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    for (int h2 = h_P_beg; h2 <= h_P_end; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      double rho_h2 = 0.0;
      if (t2z_h2 == -1) {
        rho_h2 = rho_n;
      } else {
        rho_h2 = rho_p;
      }

      int J2ab_min = std::abs(j2_h1 - j2_h2);
      int J2ab_max = j2_h1 + j2_h2;
      int T2z = t2z_p + t2z_h1 + t2z_h2;
      for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
        int J2_min = std::abs(J2ab - j2_p);
        int J2_max = J2ab + j2_p;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          int ch_Tz = TP.FindC_Tz(J2, T2z);

          int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, p);
          int y = x; // TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
          if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
            cout << "Wrong happened @ Norm_Ord::cal_val1f3_Tz 3" << endl;
            cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
            exit(0);
          }
          // double mat_t = this->vec_Tz[ch_Tz][x][y];
          double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
          double C1 = (J2 + 1.0) / (j2_p + 1.0);
          double mat = C1 * mat_t * rho_h2;

          // cout << "\t J2 = " << J2 << "\t T2z = " << T2z
          //     << "\t J2ab = " << J2ab<<"\t j2_p = "<<j2_p<<"\t";
          //   << endl;
          /*cout << "h1 = " << h1 << "\t h2 = " << h2 << "\t p = " <<
          p<<endl; cout << "C1 = " << C1 << "\t mat_t = " << mat_t <<"\t
          rho1 =
          "<<rho_h1<<"\t rho2 = "<<rho_h2<<endl; cout<< "\t mat = " <<
          mat<< endl; /**/
          res2 += mat;
        }
      }
    }
  }
  res = 0.5 * (res1 + res2);
  return res;
}

double Norm_Ord::cal2f3_Tz_C(State_DP &T2P_i, State_DP &T2P_f, int h_core) {
  double res = 0.0;
  int T2ab_z_i = T2P_i.T_2z;
  int T2ab_z_f = T2P_f.T_2z;
  int J2ab_i = T2P_i.J_2;
  int J2ab_f = T2P_f.J_2;
  if (T2ab_z_i != T2ab_z_f)
    return res;
  if (J2ab_i != J2ab_f)
    return res;
  int T2ab_z = T2ab_z_i;
  int J2ab = J2ab_i;
  int a = T2P_i.a;
  int b = T2P_i.b;
  int d = T2P_f.a;
  int e = T2P_f.b;

  // cout << "--- -- a = " << a << "\t b = " << b << "\t\t d = " << d
  //     << "\t e = " << e << endl;

  // double result = 0.0;
  for (int h = 0; h <= h_core; h++) {
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2max = j2_h + J2ab;
    int J2min = std::abs(j2_h - J2ab);
    int T2z = T2ab_z + t2z_h;
    for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) {
      int ch_Tz = TP.FindC_Tz(J2, T2z);
      /*cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " <<
      J2ab
           << endl;
      cout << "a = " << a << "\t b = " << b << "\t h = " << h <<
      endl;/**/
      /*if (SP.state_Tz[a].l == 0 && SP.state_Tz[b].l == 0 && a == b &&
      b == h) continue; if (SP.state_Tz[d].l == 0 && SP.state_Tz[e].l ==
      0 && d == e && e == h) continue;/**/
      int x = TP.FindTP_Tz(J2, T2z, J2ab, a, b, h);
      // cout << "d = " << d << "\t e = " << e << "\t h = " << h <<
      // endl;/**/
      int y = TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
      if (x < 0) {
        cout << "Wrong at cal2f3_Tz x< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
        exit(0);
      }
      if (y < 0) {
        cout << "Wrong at cal2f3_Tz y< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "d = " << d << "\t e = " << e << "\t h = " << h << endl;
        exit(0);
      }
      // cout << "ch_Tz = " << ch_Tz << "\t x = " << x << "\t y = " << y
      // << endl;
      if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
        cout << "Wrong happened @ Norm_Ord::cal_val2f3_Tz 1" << endl;
        cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
        exit(0);
      }
      // double mat_t = this->vec_Tz[ch_Tz][x][y];
      double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
      // cout << "mat_t = " << mat_t << endl;
      double mat = (J2 + 1.0) / (J2ab + 1.0) * mat_t;
      res += mat;
    }
  }
  return res;
}

double Norm_Ord::cal2f3_Tz_P(State_DP &T2P_i, State_DP &T2P_f, int h_P_beg,
                             int h_P_end, vector<double> rho_p_vec_t,
                             vector<double> rho_n_vec_t) {
  double res = 0.0;
  int T2ab_z_i = T2P_i.T_2z;
  int T2ab_z_f = T2P_f.T_2z;
  int J2ab_i = T2P_i.J_2;
  int J2ab_f = T2P_f.J_2;
  if (T2ab_z_i != T2ab_z_f)
    return res;
  if (J2ab_i != J2ab_f)
    return res;
  int T2ab_z = T2ab_z_i;
  int J2ab = J2ab_i;
  int a = T2P_i.a;
  int b = T2P_i.b;
  int d = T2P_f.a;
  int e = T2P_f.b;

  double rho_p, rho_n;

  switch (T2ab_z) {
  case -2:
    // T2z = -1; nn
    rho_p = rho_p_vec_t[0];
    rho_n = rho_n_vec_t[0];
    break;

  case 0:
    // T2z = 0; pn
    rho_p = rho_p_vec_t[1];
    rho_n = rho_n_vec_t[1];
    break;

  case 2:
    // T2z = 1; pp
    rho_p = rho_p_vec_t[2];
    rho_n = rho_n_vec_t[2];
    break;

  default:
    cout << "Wrong 0 at Wrong happen at Norm_Ord::cal2f3_Tz_P " << endl;
    exit(0);
  }

  /*cout << "--- -- a = " << a << "\t b = " << b << "\t\t d = " << d
       << "\t e = " << e << endl;/**/

  // double result = 0.0;
  for (int h = h_P_beg; h <= h_P_end; h++) {
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    double rho = 0.0;
    if (t2z_h == 1) {
      // sum over porton
      rho = rho_p;
    } else if (t2z_h == -1) {
      // sum over neutron
      rho = rho_n;
    } else {
      cout << "Wrong happen at Norm_Ord::cal2f3_Tz_P" << endl;
      exit(0);
    }

    int J2max = j2_h + J2ab;
    int J2min = std::abs(j2_h - J2ab);
    int T2z = T2ab_z + t2z_h;
    for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) {
      int ch_Tz = TP.FindC_Tz(J2, T2z);
      /*cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " <<
      J2ab
           << endl;
      //cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
      /**/
      if (SP.state_Tz[a].l == 0 && SP.state_Tz[b].l == 0 && a == b && b == h)
        continue;
      if (SP.state_Tz[d].l == 0 && SP.state_Tz[e].l == 0 && d == e && e == h)
        continue;
      int x = TP.FindTP_Tz(J2, T2z, J2ab, a, b, h);
      // cout << "d = " << d << "\t e = " << e << "\t h = " << h <<
      // endl;/**/
      int y = TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
      if (x < 0) {
        cout << "Wrong at cal2f3_Tz x< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
        exit(0);
      }
      if (y < 0) {
        cout << "Wrong at cal2f3_Tz y< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "d = " << d << "\t e = " << e << "\t h = " << h << endl;
        exit(0);
      }
      // cout << "ch_Tz = " << ch_Tz << "\t x = " << x << "\t y = " << y
      // << endl;
      if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
        cout << "Wrong happened @ Norm_Ord::cal_val2f3_Tz  2" << endl;
        cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
        exit(0);
      }
      // double mat_t = this->vec_Tz[ch_Tz][x][y];
      double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
      // cout << "mat_t = " << mat_t << endl;
      double mat = (J2 + 1.0) / (J2ab + 1.0) * mat_t * rho;
      // cout << "h = "<<h<<"\t mat_t = " << mat_t <<"\t rho = "<<rho
      // <<"\t mat = " << mat << endl;

      res += mat;
    }
  }
  // cout<<"\t !!!  result = "<<res<<endl;
  return res;
}

double Norm_Ord::cal2f3_Tz_P2(State_DP &T2P_i, State_DP &T2P_f, int h_P_beg,
                              int h_P_end, vector<vector<double>> rho_val) {
  double res = 0.0;
  int T2ab_z_i = T2P_i.T_2z;
  int T2ab_z_f = T2P_f.T_2z;
  int J2ab_i = T2P_i.J_2;
  int J2ab_f = T2P_f.J_2;
  if (T2ab_z_i != T2ab_z_f)
    return res;
  if (J2ab_i != J2ab_f)
    return res;
  int T2ab_z = T2ab_z_i;
  int J2ab = J2ab_i;
  int a = T2P_i.a;
  int b = T2P_i.b;
  int d = T2P_f.a;
  int e = T2P_f.b;

  double rho_p, rho_n;

  /*cout << "--- -- a = " << a << "\t b = " << b << "\t\t d = " << d
       << "\t e = " << e << endl;/**/

  // double result = 0.0;
  for (int h = h_P_beg; h <= h_P_end; h++) {
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int Tz_index = T2ab_z / 2 + 1;
    double rho = rho_val[Tz_index][h];

    int J2max = j2_h + J2ab;
    int J2min = std::abs(j2_h - J2ab);
    int T2z = T2ab_z + t2z_h;
    for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) {
      int ch_Tz = TP.FindC_Tz(J2, T2z);
      /*cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " <<
      J2ab
           << endl;
      //cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
      /**/
      if (SP.state_Tz[a].l == 0 && SP.state_Tz[b].l == 0 && a == b && b == h)
        continue;
      if (SP.state_Tz[d].l == 0 && SP.state_Tz[e].l == 0 && d == e && e == h)
        continue;
      int x = TP.FindTP_Tz(J2, T2z, J2ab, a, b, h);
      // cout << "d = " << d << "\t e = " << e << "\t h = " << h <<
      // endl;/**/
      int y = TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
      if (x < 0) {
        cout << "Wrong at cal2f3_Tz x< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
        exit(0);
      }
      if (y < 0) {
        cout << "Wrong at cal2f3_Tz y< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "d = " << d << "\t e = " << e << "\t h = " << h << endl;
        exit(0);
      }
      // cout << "ch_Tz = " << ch_Tz << "\t x = " << x << "\t y = " << y
      // << endl;
      if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
        cout << "Wrong happened @ Norm_Ord::cal_val2f3_Tz 3" << endl;
        cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
        cout << "ch_Tz : " << ch_Tz << "\t x : " << x << "\t y : " << y << endl;
        cout << "a : " << a << "\t b : " << b << "\t h : " << h
             << "\t d : " << d << "\t e : " << e << endl;
        exit(0);
      }
      // double mat_t = this->vec_Tz[ch_Tz][x][y];
      double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
      // cout << "mat_t = " << mat_t << endl;
      double mat = (J2 + 1.0) / (J2ab + 1.0) * mat_t * rho;
      /*if (T2ab_z == 2) {
        cout << "T2ab_z = " << T2ab_z << "\t ";
        cout << "h = " << h << "\t mat_t = " << mat_t << "\t rho = " <<
      rho
             << "\t mat = " << mat << endl;
      }/**/
      res += mat;
    }
  }
  // cout<<"\t !!!  result = "<<res<<endl;
  return res;
}

void Norm_Ord::cal_val0f3_Tz(int h_core, int h_P_beg, int h_P_end, int Num_p,
                             int Num_n) {

  double res = 0.0;
  int Num_C_p = 0;
  int Num_C_n = 0;

  for (int h = 0; h <= h_core; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_C_n += j2 + 1;
    } else if (t2z == 1) {
      Num_C_p += j2 + 1;
    } else {
      cout << "Wrong 1 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "Particles Number, Proton : " << Num_p << "\t Neutron : " << Num_n
       << endl;
  cout << "Particles Number in core, Proton : " << Num_C_p
       << "\t Neutron : " << Num_C_n << endl;

  int Num_Full_P_p = 0;
  int Num_Full_P_n = 0;

  for (int h = h_P_beg; h <= h_P_end; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_Full_P_n += j2 + 1;
    } else if (t2z == 1) {
      Num_Full_P_p += j2 + 1;
    } else {
      cout << "Wrong 2 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "obital 2j+1 , Proton : " << Num_Full_P_p
       << "\t Neutron : " << Num_Full_P_n << endl;

  rho_p_H0 = (1.0 * (Num_p - Num_C_p)) / (1.0 * Num_Full_P_p);
  rho_n_H0 = (1.0 * (Num_n - Num_C_n)) / (1.0 * Num_Full_P_n);

  for (int h1 = h_P_beg; h1 <= h_P_end; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    double rho_h1 = 0.0;
    if (t2z_h1 == -1) {
      rho_h1 = rho_n_H0;
    } else {
      rho_h1 = rho_p_H0;
    }
    if (rho_h1 == 0.0)
      continue;

    for (int h2 = h_P_beg; h2 <= h_P_end; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      double rho_h2 = 0.0;
      if (t2z_h2 == -1) {
        rho_h2 = rho_n_H0;
      } else {
        rho_h2 = rho_p_H0;
      }
      if (rho_h2 == 0.0)
        continue;

      for (int h3 = h_P_beg; h3 <= h_P_end; h3++) {
        int j2_h3 = SP.state_Tz[h3].j_2;
        int t2z_h3 = SP.state_Tz[h3].tz_2;
        double rho_h3 = 0.0;
        if (t2z_h3 == -1) {
          rho_h3 = rho_n_H0;
        } else {
          rho_h3 = rho_p_H0;
        }
        if (rho_h3 == 0.0)
          continue;

        int J2ab_min = std::abs(j2_h1 - j2_h2);
        int J2ab_max = j2_h1 + j2_h2;
        int T2z = t2z_h1 + t2z_h2 + t2z_h3;
        for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
          int J2_min = std::abs(J2ab - j2_h3);
          int J2_max = J2ab + j2_h3;
          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            int ch_Tz = TP.FindC_Tz(J2, T2z);

            int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, h3);
            int y = x; // TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
            if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
              cout << "Wrong happened @ Norm_Ord::cal_val0f3_Tz" << endl;
              cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
              exit(0);
            }
            // double mat_t = this->.val[ch_Tz][x][y];
            double mat_t = this->vec_Tz_S.val(ch_Tz, x, y);
            double C1 = (J2 + 1.0);
            double mat = C1 * mat_t * rho_h1 * rho_h2 * rho_h3;
            res += mat;
          }
        }
      }
    }
  }
  val_0f3 = res * 1.0 / 6.0;
}

void Norm_Ord::build_vec1f3_Tz_C(int h_core, int p_beg, int p_end) {

  int size = SP.state_Tz.size();
  vec_1f3_C.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3_C[i] = 0.0;
  }
  for (int i = p_beg; i <= p_end; i++) {
    // cout << "---- i = " << i << endl;
    vec_1f3_C[i] = this->cal1f3_Tz_C(SP.state_Tz[i], h_core);
  }
}

void Norm_Ord::build_vec1f3_Tz_P(int p_beg, int p_end, int h_P_beg, int h_P_end,
                                 vector<double> rho_p_sp_t,
                                 vector<double> rho_n_sp_t) {
  int size = SP.state_Tz.size();
  vec_1f3_P.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3_P[i] = 0.0;
  }
  for (int i = p_beg; i <= p_end; i++) {
    // cout << "---- i = " << i << endl;
    vec_1f3_P[i] = this->cal1f3_Tz_P(SP.state_Tz[i], h_P_beg, h_P_end,
                                     rho_p_sp_t, rho_n_sp_t);
  }
}
void Norm_Ord::build_vec1f3_Tz_1C1P(int h_core, int p_beg, int p_end,
                                    int h_P_beg, int h_P_end,
                                    vector<double> rho_p_sp_t,
                                    vector<double> rho_n_sp_t) {
  int size = SP.state_Tz.size();
  vec_1f3_1C1P.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3_1C1P[i] = 0.0;
  }
  for (int i = p_beg; i <= p_end; i++) {
    // cout << "---- i = " << i << endl;
    vec_1f3_1C1P[i] = this->cal1f3_Tz_1C1P(SP.state_Tz[i], h_core, h_P_beg,
                                           h_P_end, rho_p_sp_t, rho_n_sp_t);
  }
}

void Norm_Ord::build_vec1f3_Tz(int h_core, int p_beg, int p_end, int h_P_beg,
                               int h_P_end, int Num_p, int Num_n) {
  int Num_C_p = 0;
  int Num_C_n = 0;

  for (int h = 0; h <= h_core; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_C_n += j2 + 1;
    } else if (t2z == 1) {
      Num_C_p += j2 + 1;
    } else {
      cout << "Wrong 1 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "Particles Number, Proton : " << Num_p << "\t Neutron : " << Num_n
       << endl;
  cout << "Particles Number in core, Proton : " << Num_C_p
       << "\t Neutron : " << Num_C_n << endl;

  int Num_Full_P_p = 0;
  int Num_Full_P_n = 0;

  for (int h = h_P_beg; h <= h_P_end; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_Full_P_n += j2 + 1;
    } else if (t2z == 1) {
      Num_Full_P_p += j2 + 1;
    } else {
      cout << "Wrong 2 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "obital 2j+1 , Proton : " << Num_Full_P_p
       << "\t Neutron : " << Num_Full_P_n << endl;
  this->Num_valence_p = Num_p - Num_C_p;
  this->Num_valence_n = Num_n - Num_C_n;
  cout << "\t Num_valence_p = " << Num_valence_p
       << "\t Num_valence_n = " << Num_valence_n << endl;

  rho_p_vec_sp.resize(2); // tz = -1, 1.  n, p
  rho_n_vec_sp.resize(2); // tz = -1, 1.  n, p

  for (int i = 0; i < 2; i++) {
    int active_p;
    int active_n;
    switch (i) {
    case 0:
      // tz = -1; n
      active_n = 1;
      active_p = 0;
      break;

    case 1:
      // tz = 0; pn
      active_n = 0;
      active_p = 1;
      break;
    }

    double rho_p = (1.0 * (Num_p - Num_C_p - active_p)) / (1.0 * Num_Full_P_p);
    double rho_n = (1.0 * (Num_n - Num_C_n - active_n)) / (1.0 * Num_Full_P_n);

    if (rho_p < 0)
      rho_p = 0.0;
    if (rho_n < 0)
      rho_n = 0.0;

    rho_p_vec_sp[i] = rho_p;
    rho_n_vec_sp[i] = rho_n;
  }

  cout << "Proton rho_sp tz=-1 : " << rho_p_vec_sp[0]
       << "\t tz=1 : " << rho_p_vec_sp[1] << endl;
  cout << "Neutron rho_sp tz=-1 : " << rho_n_vec[0]
       << "\t tz=1 : " << rho_n_vec_sp[1] << endl;

  this->build_vec1f3_Tz_C(h_core, p_beg, p_end);

  this->build_vec1f3_Tz_P(p_beg, p_end, h_P_beg, h_P_end, rho_p_vec_sp,
                          rho_n_vec_sp);

  this->build_vec1f3_Tz_1C1P(h_core, p_beg, p_end, h_P_beg, h_P_end,
                             rho_p_vec_sp, rho_n_vec_sp);

  int size = SP.state_Tz.size();
  vec_1f3.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3[i] = vec_1f3_C[i] + vec_1f3_P[i] + vec_1f3_1C1P[i];
  }

  /**/
}

void Norm_Ord::build_vec1f3_Tz2(int h_core, int p_beg, int p_end, int h_P_beg,
                                int h_P_end, int Num_p, int Num_n, int flag_P) {
  int Num_C_p = 0;
  int Num_C_n = 0;

  int size = SP.state_Tz.size();

  vec_1f3_P.resize(size);
  vec_1f3_1C1P.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3_P[i] = 0.0;
    vec_1f3_1C1P[i] = 0.0;
  }

  this->build_vec1f3_Tz_C(h_core, p_beg, p_end);

  vec_1f3.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3[i] = vec_1f3_C[i] + vec_1f3_P[i] + vec_1f3_1C1P[i];
  }

  /**/
}
void Norm_Ord::build_vec1f3_Tz3(int sp_core_Tz, int sp_P_beg_Tz,
                                int sp_P_end_Tz) {
  int size = SP.state_Tz.size();
  vec_1f3_P.resize(size);
  vec_1f3_1C1P.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3_P[i] = 0.0;
    vec_1f3_1C1P[i] = 0.0;
  }

  this->build_vec1f3_Tz_C(sp_core_Tz, sp_P_beg_Tz, sp_P_end_Tz);
  vec_1f3.resize(size);
  for (int i = 0; i < size; i++) {
    vec_1f3[i] = vec_1f3_C[i];
  }
}

void Norm_Ord::build_vec2f3_Tz_C(int h_core) {

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3_C.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3_C[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3_C[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3_C[ch][i][j] = 0.0;
      }
    }
  }
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double mat =
            this->cal2f3_Tz_C(state_T2P_Tz[ch][i], state_T2P_Tz[ch][j], h_core);
        vec_2f3_C[ch][i][j] = mat;
      }
    }
  }
}

void Norm_Ord::build_vec2f3_Tz_P(int h_P_beg, int h_P_end,
                                 vector<double> rho_p_vec_t,
                                 vector<double> rho_n_vec_t) {

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3_P.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3_P[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3_P[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3_P[ch][i][j] = 0.0;
      }
    }
  }
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double mat =
            this->cal2f3_Tz_P(state_T2P_Tz[ch][i], state_T2P_Tz[ch][j], h_P_beg,
                              h_P_end, rho_p_vec_t, rho_n_vec_t);
        vec_2f3_P[ch][i][j] = mat;
      }
    }
  } /**/ /*
   for (int ch = 0; ch < size_ch_Tz; ch++) {
     int size = state_T2P_Tz[ch].size();
     for (int i = 0; i < 1; i++) {
       for (int j = 0; j < 1; j++) {
         double mat =
             this->cal2f3_Tz_P(state_T2P_Tz[ch][i], state_T2P_Tz[ch][j],
   h_P_beg,
                               h_P_end, rho_p_vec_t, rho_n_vec_t);
         // vec_2f3_P[ch][i][j] = mat;
       }
     }
   }/**/
}

void Norm_Ord::build_vec2f3_Tz_P2(int h_P_beg, int h_P_end,
                                  vector<vector<double>> rho_val) {

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3_P.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3_P[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3_P[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3_P[ch][i][j] = 0.0;
      }
    }
  }
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double mat =
            this->cal2f3_Tz_P2(state_T2P_Tz[ch][i], state_T2P_Tz[ch][j],
                               h_P_beg, h_P_end, rho_val);
        vec_2f3_P[ch][i][j] = mat;
      }
    }
  } /**/ /*
   for (int ch = 0; ch < size_ch_Tz; ch++) {
     int size = state_T2P_Tz[ch].size();
     for (int i = 0; i < 1; i++) {
       for (int j = 0; j < 1; j++) {
         double mat =
             this->cal2f3_Tz_P(state_T2P_Tz[ch][i], state_T2P_Tz[ch][j],
   h_P_beg,
                               h_P_end, rho_p_vec_t, rho_n_vec_t);
         // vec_2f3_P[ch][i][j] = mat;
       }
     }
   }/**/
}

void Norm_Ord::build_vec2f3_Tz(int h_core, int h_P_beg, int h_P_end, int Num_p,
                               int Num_n) {
  int Num_C_p = 0;
  int Num_C_n = 0;

  for (int h = 0; h <= h_core; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_C_n += j2 + 1;
    } else if (t2z == 1) {
      Num_C_p += j2 + 1;
    } else {
      cout << "Wrong 1 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "Particles Number, Proton : " << Num_p << "\t Neutron : " << Num_n
       << endl;
  cout << "Particles Number in core, Proton : " << Num_C_p
       << "\t Neutron : " << Num_C_n << endl;

  int Num_Full_P_p = 0;
  int Num_Full_P_n = 0;

  for (int h = h_P_beg; h <= h_P_end; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_Full_P_n += j2 + 1;
    } else if (t2z == 1) {
      Num_Full_P_p += j2 + 1;
    } else {
      cout << "Wrong 2 happened at build_vec2f3_Tz" << endl;
      exit(0);
    }
  }
  cout << "obital 2j+1 , Proton : " << Num_Full_P_p
       << "\t Neutron : " << Num_Full_P_n << endl;

  rho_p_vec.resize(3); // tz = -1, 0, 1.  nn, pn, pp
  rho_n_vec.resize(3); // tz = -1, 0, 1.  nn, pn, pp

  for (int i = 0; i < 3; i++) {
    int active_p;
    int active_n;
    switch (i) {
    case 0:
      // tz = -1; nn
      active_n = 2;
      active_p = 0;
      break;

    case 1:
      // tz = 0; pn
      active_n = 1;
      active_p = 1;
      break;

    case 2:
      // tz = 1; pp
      active_n = 0;
      active_p = 2;
      break;
    }

    double rho_p = (1.0 * (Num_p - Num_C_p - active_p)) / (1.0 * Num_Full_P_p);
    double rho_n = (1.0 * (Num_n - Num_C_n - active_n)) / (1.0 * Num_Full_P_n);

    if (rho_p < 0)
      rho_p = 0.0;
    if (rho_n < 0)
      rho_n = 0.0;

    rho_p_vec[i] = rho_p;
    rho_n_vec[i] = rho_n;
  }

  cout << "Proton rho tz=-1 : " << rho_p_vec[0] << "\t tz=0 : " << rho_p_vec[1]
       << "\t tz=1 : " << rho_p_vec[2] << endl;
  cout << "Neutron rho tz=-1 : " << rho_n_vec[0] << "\t tz=0 : " << rho_n_vec[1]
       << "\t tz=1 : " << rho_n_vec[2] << endl;

  this->build_vec2f3_Tz_P(h_P_beg, h_P_end, rho_p_vec, rho_n_vec);

  this->build_vec2f3_Tz_C(h_core);

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3[ch][i][j] = vec_2f3_C[ch][i][j] + vec_2f3_P[ch][i][j];
      }
    }
  } /**/
}

void Norm_Ord::build_vec2f3_Tz2(int h_core, int h_P_beg, int h_P_end, int Num_p,
                                int Num_n, int flag_P) {
  int Num_C_p = 0;
  int Num_C_n = 0;

  for (int h = 0; h <= h_core; h++) {
    int j2 = SP.state_Tz[h].j_2;
    int t2z = SP.state_Tz[h].tz_2;
    if (t2z == -1) {
      Num_C_n += j2 + 1;
    } else if (t2z == 1) {
      Num_C_p += j2 + 1;
    } else {
      cout << "Wrong 1 happened at build_vec2f3_Tz2" << endl;
      exit(0);
    }
  }
  cout << "Particles Number, Proton : " << Num_p << "\t Neutron : " << Num_n
       << endl;
  cout << "Particles Number in core, Proton : " << Num_C_p
       << "\t Neutron : " << Num_C_n << endl;

  int size_sp_tz = SP.state_Tz.size();

  int nuclei_num_max = 6;                            // For p-shell
  vector<vector<vector<double>>> num_rho_ingredient; // [neutron][proton][sp_tz]
  num_rho_ingredient.resize(nuclei_num_max + 1);
  for (int i = 0; i <= nuclei_num_max; i++) {
    num_rho_ingredient[i].resize(nuclei_num_max + 1);
    for (int j = 0; j <= nuclei_num_max; j++) {
      num_rho_ingredient[i][j].resize(size_sp_tz);
      for (int k = 0; k < size_sp_tz; k++) {
        num_rho_ingredient[i][j][k] = -1;
      }
    }
  }
  // 8Be density, 2 neutron 2 proton in valence space
  for (int i = h_P_beg; i <= h_P_end; i++) {
    int tz_2 = SP.state_Tz[i].tz_2;
    int j_2 = SP.state_Tz[i].j_2;
    int l = SP.state_Tz[i].l;
    if (tz_2 == 1) {
      // proton
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[2][2][i] = 1.4;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[2][2][i] = 0.6;
      }
    }
    if (tz_2 == -1) {
      // neutron
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[2][2][i] = 1.4;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[2][2][i] = 0.6;
      }
    }
  }
  // 8B density, 1 neutron 3 proton in valence space
  for (int i = h_P_beg; i <= h_P_end; i++) {
    int tz_2 = SP.state_Tz[i].tz_2;
    int j_2 = SP.state_Tz[i].j_2;
    int l = SP.state_Tz[i].l;
    if (tz_2 == 1) {
      // proton
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[1][3][i] = 2.8;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[1][3][i] = 0.2;
      }
    }
    if (tz_2 == -1) {
      // neutron
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[1][3][i] = 0.9;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[1][3][i] = 0.1;
      }
    }
  }
  // 8Li density, 3 neutron 1 proton in valence space
  for (int i = h_P_beg; i <= h_P_end; i++) {
    int tz_2 = SP.state_Tz[i].tz_2;
    int j_2 = SP.state_Tz[i].j_2;
    int l = SP.state_Tz[i].l;
    if (tz_2 == 1) {
      // proton
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[3][1][i] = 0.9;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[3][1][i] = 0.1;
      }
    }
    if (tz_2 == -1) {
      // neutron
      if (l == 1 && j_2 == 3) {
        num_rho_ingredient[3][1][i] = 2.8;
      }
      if (l == 1 && j_2 == 1) {
        num_rho_ingredient[3][1][i] = 0.2;
      }
    }
  }

  this->rho_vec.resize(3);

  for (int i = 0; i < 3; i++) {
    rho_vec[i].resize(size_sp_tz);
    for (int j = 0; j < size_sp_tz; j++) {
      rho_vec[i][j] = 0.0;
    }
  }

  for (int i = 0; i < 3; i++) {
    int active_p;
    int active_n;
    switch (i) {
    case 0:
      // tz = -1; nn
      active_n = 2;
      active_p = 0;
      break;

    case 1:
      // tz = 0; pn
      active_n = 1;
      active_p = 1;
      break;

    case 2:
      // tz = 1; pp
      active_n = 0;
      active_p = 2;
      break;
    }

    int num_valence_p = (Num_p - Num_C_p - active_p);
    int num_valence_n = (Num_n - Num_C_n - active_n);

    for (int j = h_P_beg; j <= h_P_end; j++) {
      int tz_2 = SP.state_Tz[j].tz_2;
      int j_2 = SP.state_Tz[j].j_2;
      int Num_Orb = j_2 + 1;
      double num_rho = num_rho_ingredient[num_valence_n][num_valence_p][j];
      if (num_rho < 0) {
        cout << "Wrong at build_vec2f3_Tz2 num_rho<0" << endl;
        cout << "j = " << j << endl;
        exit(0);
      }
      double rho = num_rho / (1.0 * Num_Orb);
      rho_vec[i][j] = rho;
    }
  }

  for (int i = 0; i < 3; i++) {
    int Tz = i - 1;
    for (int j = 0; j < size_sp_tz; j++) {
      int j_2 = SP.state_Tz[j].j_2;
      int tz_2 = SP.state_Tz[j].tz_2;
      if (rho_vec[i][j] > 0) {
        cout << "Tz : " << Tz << "\t\t";
        cout << "j_2 = " << j_2 << "\t ";
        cout << "tz_2 = " << tz_2 << "\t ";
        cout << "rho = " << rho_vec[i][j] << endl;
      }
    }
  }
  if (flag_P == 1) {
    this->build_vec2f3_Tz_P2(h_P_beg, h_P_end, rho_vec);
  } else {
    int size_ch_Tz = channel_T2P_Tz.size();
    vec_2f3_P.resize(size_ch_Tz);
    for (int ch = 0; ch < size_ch_Tz; ch++) {
      int size = state_T2P_Tz[ch].size();
      vec_2f3_P[ch].resize(size);
      for (int i = 0; i < size; i++) {
        vec_2f3_P[ch][i].resize(size);
        for (int j = 0; j < size; j++) {
          vec_2f3_P[ch][i][j] = 0;
        }
      }
    } /**/
  }

  this->build_vec2f3_Tz_C(h_core);

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3[ch][i][j] = vec_2f3_C[ch][i][j] + vec_2f3_P[ch][i][j];
      }
    }
  } /**/
} /**/

void Norm_Ord::build_vec2f3_Tz3(int sp_core_Tz) {
  this->build_vec2f3_Tz_C(sp_core_Tz);

  int size_ch_Tz = channel_T2P_Tz.size();
  vec_2f3.resize(size_ch_Tz);
  for (int ch = 0; ch < size_ch_Tz; ch++) {
    int size = state_T2P_Tz[ch].size();
    vec_2f3[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2f3[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2f3[ch][i][j] = vec_2f3_C[ch][i][j];
      }
    }
  } /**/
}

void Norm_Ord::print_T2P_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "a b JJ-coupling state_Tz" << endl;
  // file << "Energy cut-off is : " << Emax << " ( hbar-omega )" <<
  // endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state_T2P_Tz.size(); i++) {
    total_num += state_T2P_Tz[i].size();
  }
  file << total_num << endl;
  file << channel_T2P_Tz.size() << "\t channel_Tz number" << endl;
  for (int i = 0; i < channel_T2P_Tz.size(); i++) {
    file << state_T2P_Tz[i].size() << "\t";
  }
  file << endl;
  for (int i = 0; i < channel_T2P_Tz.size(); i++) {
    file << "===== channel_Tz : =====" << i << endl;
    file << "T2ab_z \t"
         << "J2ab \t"
         << "a \t"
         << "b \t"
         << "E \t"
         << "index" << endl;
    for (int j = 0; j < state_T2P_Tz[i].size(); j++) {
      file << state_T2P_Tz[i][j].T_2z << "\t";
      file << state_T2P_Tz[i][j].J_2 << "\t";
      file << state_T2P_Tz[i][j].a << "\t";
      file << state_T2P_Tz[i][j].b << "\t";
      file << state_T2P_Tz[i][j].E << "\t";
      file << state_T2P_Tz[i][j].index << "\t";
      file << endl;
    }
  }
  file.close();
}

void Norm_Ord::print_channel_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel_T2P_Tz.size();
  file << " TP channel_Tz Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J_2 \t"
       << "T_2z \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel_T2P_Tz[i].num << "\t" << channel_T2P_Tz[i].J_2
         << "\t" << channel_T2P_Tz[i].T_2z << "\t" << channel_T2P_Tz[i].E_min
         << "\t" << channel_T2P_Tz[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel_T2P_Tz[i].indexE_beg.size(); j++) {
      file << channel_T2P_Tz[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel_T2P_Tz[i].indexE.size(); j++) {
      file << channel_T2P_Tz[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}

void Norm_Ord::print_vec2f3_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
          file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
}
void Norm_Ord::print_vec1f3_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = vec_1f3_C.size();
  for (int i = 0; i < size; i++) {

    file << i << "\t" << vec_1f3_C[i] << endl;
  }

  file.close();
}

void Norm_Ord::print_norm_ME_CP(std::string filename, int h_core_Tz,
                                int p_beg_Tz, int p_end_Tz) {
  std::ofstream file;
  file.open(filename);
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t ";
  for (int i = p_end; i >= p_beg; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering zero-body Energy === " << endl;
  file << "\tzero-body H0 enery : " << val_0f3 << " (MeV)" << endl;
  file << endl;

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "\tdensity in p3/2: " << endl;
  file << "\t MT\t Proton\t Neutron" << endl;
  for (int i = 0; i < 2; i++) {
    file << "\t" << i * 2 - 1 << "\t" << rho_p_vec_sp[i] << "\t"
         << rho_n_vec_sp[i] << endl;
  }
  file << "index\t n\t l\t 2*j\t 2*tz\t  E{core}   +   E{valence}  +  "
          "E{(core)(valence)}     =     E(MeV)"
       << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << std::floor(i / 2.0) + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_C[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_P[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_1C1P[i] << "\t\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3[i] << "\t ";
    file << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core + p3/2} "
          "orbits=== "
       << endl;
  file << "\tdensity in p3/2: " << endl;
  file << "\t MT\t Proton\t Neutron" << endl;
  for (int i = 0; i < 3; i++) {
    file << "\t" << i - 1 << "\t" << rho_p_vec[i] << "\t" << rho_n_vec[i]
         << endl;
  }
  file << "\t MT=-1 : <nn|V|nn>; MT=1 <pp|V|pp>; MT=0: <pn|V|pn>" << endl;
  file << "\t <(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME{core}  +  ME{valence}    =    "
          "ME (MeV)"
       << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        if (T2ab_z == 0 && t2a_p < 0)
          continue;
        if (T2ab_z == 0 && t2a < 0)
          continue; /* Only considering <pn|V|pn> */
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b + 1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t ";
        file << std::setw(8) << std::setprecision(5) << vec_2f3_C[ch][i][j]
             << "\t ";
        file << std::setw(8) << std::setprecision(5) << vec_2f3_P[ch][i][j]
             << "\t "; /**/
        file << std::setw(8) << std::setprecision(5) << vec_2f3[ch][i][j]
             << endl;

        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::print_norm_ME_CP2(std::string filename, int h_core_Tz,
                                 int p_beg_Tz, int p_end_Tz) {
  std::ofstream file;
  file.open(filename);
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t ";
  for (int i = p_end; i >= p_beg; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering zero-body Energy === " << endl;
  file << "\tzero-body H0 enery : " << 0.0 << " (MeV)" << endl;
  file << endl;

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "index\t n\t l\t 2*j\t 2*tz\t  E{core}   +   E{valence}  +  "
          "E{(core)(valence)}     =     E(MeV)"
       << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << std::floor(i / 2.0) + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_C[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_P[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_1C1P[i] << "\t\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3[i] << "\t ";
    file << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core + p3/2} "
          "orbits=== "
       << endl;
  file << "\tdensity in p3/2 & p1/2: " << endl;
  file << "\t MTab\t j_2\t tz_2\t\t rho" << endl;
  for (int i = 0; i < 3; i++) {
    int Tz = i - 1;
    for (int j = 0; j < rho_vec[i].size(); j++) {
      int j_2 = SP.state_Tz[j].j_2;
      int tz_2 = SP.state_Tz[j].tz_2;
      if (rho_vec[i][j] > 0) {
        file << "\t " << Tz << "\t ";
        file << j_2 << "\t ";
        file << tz_2 << "\t\t ";
        file << rho_vec[i][j] << endl;
      }
    }
  }
  file << "\t MT=-1 : <nn|V|nn>; MT=1 <pp|V|pp>; MT=0: <pn|V|pn>" << endl;
  file << "\t <(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME{core}  +  ME{valence}    =    "
          "ME (MeV)"
       << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        if (T2ab_z == 0 && t2a_p < 0)
          continue;
        if (T2ab_z == 0 && t2a < 0)
          continue; /* Only considering <pn|V|pn> */
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b + 1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t ";
        file << std::setw(8) << std::setprecision(5) << vec_2f3_C[ch][i][j]
             << "\t ";
        file << std::setw(8) << std::setprecision(5) << vec_2f3_P[ch][i][j]
             << "\t "; /**/
        file << std::setw(8) << std::setprecision(5) << vec_2f3[ch][i][j]
             << endl;

        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::print_with_SecOrd(std::string filename, int h_core_Tz,
                                 int p_beg_Tz, int p_end_Tz,
                                 vector<double> vec_sec,
                                 vector<vector<double>> vec_sec_sub, double hw,
                                 int reg_pow, int Lambda_Cut, int E_cut) {
  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t cE = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hw << endl;
  file << "\t The Power of Regulator : " << reg_pow << endl;
  file << " \t +++++++++++++++++" << endl;
  file << "\t    E_cut = " << E_cut << endl;
  file << " \t +++++++++++++++++" << endl;
  file << " ----------------" << endl;

  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t ";
  for (int i = p_end; i >= p_beg; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering zero-body Energy === " << endl;
  file << "\tzero-body H0 enery : " << 0.0 << " (MeV)" << endl;
  file << endl;

  file << " === 13 Second Order Diagrams value === " << endl;
  file << "a\t A1\t\t A2\t\t B1\t\t B2\t\t C1\t\t C2\t\t D1\t\t D2\t\t E1\t\t "
          "E2\t\t E3\t\t F1\t\t "
          "F2\t\t  "
       << endl;
  cout << "\t vec_sec_sub.size : " << vec_sec_sub.size() << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    /*file << std::floor(i / 2.0) + 1 << "\t";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[0][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[1][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[2][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[3][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[4][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[5][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[6][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[7][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[8][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[9][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[10][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[11][i] << "\t ";
    file << std::setw(6) << std::setprecision(4) << vec_sec_sub[12][i] << "\t "
         << endl;/**/
    file << std::floor(i / 2.0) + 1 << "\t";
    file << std::setw(6) << vec_sec_sub[0][i] << "\t";
    file << std::setw(6) << vec_sec_sub[1][i] << "\t";
    file << std::setw(6) << vec_sec_sub[2][i] << "\t";
    file << std::setw(6) << vec_sec_sub[3][i] << "\t";
    file << std::setw(6) << vec_sec_sub[4][i] << "\t";
    file << std::setw(6) << vec_sec_sub[5][i] << "\t";
    file << std::setw(6) << vec_sec_sub[6][i] << "\t";
    file << std::setw(6) << vec_sec_sub[7][i] << "\t";
    file << std::setw(6) << vec_sec_sub[8][i] << "\t";
    file << std::setw(6) << vec_sec_sub[9][i] << "\t";
    file << std::setw(6) << vec_sec_sub[10][i] << "\t";
    file << std::setw(6) << vec_sec_sub[11][i] << "\t";
    file << std::setw(6) << vec_sec_sub[12][i] << "\t" << endl;
  }
  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "index\t n\t l\t 2*j\t 2*tz\t  E{core}   +   E{valence}  +  "
          "E{(core)(valence)}   +   Second_Order    =     E(MeV)"
       << endl;

  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << std::floor(i / 2.0) + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_C[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_P[i] << "\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3_1C1P[i] << "\t\t ";
    file << std::setw(8) << std::setprecision(5) << vec_sec[i] << "\t\t ";
    file << std::setw(8) << std::setprecision(5) << vec_1f3[i] + vec_sec[i]
         << "\t ";
    file << endl;
  }
  file << endl;

  file << "\t MT=-1 : <nn|V|nn>; MT=1 <pp|V|pp>; MT=0: <pn|V|pn>" << endl;
  file << "\t <(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME{core}  " << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        if (T2ab_z == 0 && t2a_p < 0)
          continue;
        if (T2ab_z == 0 && t2a < 0)
          continue; /* Only considering <pn|V|pn> */
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b + 1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t ";
        file << std::setw(8) << std::setprecision(5) << vec_2f3_C[ch][i][j]
             << endl;

        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::print_norm_ME_C(std::string filename, int h_core_Tz,
                               int p_beg_Tz, int p_end_Tz, double hw, int reg_pow,
                               int Lambda_Cut, vector<int> oc_o2z) {
  // y= oc_o2z[x] : y is state_Tz;  x is from others
  vector<int> oc_z2o;
  int size_oc = oc_o2z.size();
  oc_z2o.resize(size_oc);
  for (int x = 0; x < size_oc; x++) {
    int y = oc_o2z[x];
    oc_z2o[y] = x;
  }

  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t gA = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hw << endl;
  file << "\t The Power of Regulator : " << reg_pow << endl;
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t";
  for (int i = p_end; i >= p_beg; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "\t Num: " << p_end_Tz - p_beg_Tz + 1 << endl;
  file << "\tindex\t n\t l\t 2*j\t 2*tz\t E (MeV)" << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    int y = std::floor(i / 2.0);
    int x = oc_z2o[y];
    file << "\t" << x + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(6) << vec_1f3_C[i] << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core} orbits=== "
       << endl;
  file << "\tMT=-1 : <nn|V|nn>; MT=1 <pp|V|pp>; MT=0: <pn|V|pn>" << endl;
  file << "\t<(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME (MeV)" << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        if (T2ab_z == 0 && t2a_p < 0)
          continue;
        if (T2ab_z == 0 && t2a < 0)
          continue; /* Only considering <pn|V|pn> */
        /*file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b + 1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";/**/
        file << oc_z2o[ap] + 1 << "\t" << oc_z2o[bp] + 1 << "\t"
             << oc_z2o[a] + 1 << "\t" << oc_z2o[b] + 1 << "\t" << J2ab / 2
             << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3_C[ch][i][j] << endl;
        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::print_norm_ME_CP_Tz(std::string filename, int h_core_Tz,
                                   int p_beg_Tz, int p_end_Tz) {
  std::ofstream file;
  file.open(filename);
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\t index\t n\t l\t 2*j\t 2*tz";
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core_Tz; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "\t index\t n\t l\t 2*j\t 2*tz\t E (MeV)" << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    // file << "\t" << std::floor(i / 2.0) + 1 << "\t";
    file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(6) << vec_1f3_C[i] << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core + p3/2} "
          "orbits=== "
       << endl;
  file << " <(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << "\t MT=-1 : nn; MT=1 pp; MT=0: pn or np" << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME (MeV)" << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        // if (T2ab_z == 0 && ap == bp && ap_Tz > 0)
        //  continue;
        // if (T2ab_z == 0 && a == b && a_Tz > 0)
        //  continue;
        file << ap_Tz + 1 << "\t" << bp_Tz + 1 << "\t" << a_Tz + 1 << "\t"
             << b_Tz + 1 << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j] << endl;
        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::print_norm_ME_C_Tz(std::string filename, int h_core_Tz,
                                  int p_beg_Tz, int p_end_Tz) {
  std::ofstream file;
  file.open(filename);
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\t index\t n\t l\t 2*j\t 2*tz";
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core_Tz; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "\t index\t n\t l\t 2*j\t 2*tz\t E (MeV)" << endl;
  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    // file << "\t" << std::floor(i / 2.0) + 1 << "\t";
    file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    file << std::setw(6) << vec_1f3_C[i] << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core} orbits=== "
       << endl;
  file << " <(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << "\t MT=-1 : nn; MT=1 pp; MT=0: pn or np" << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME (MeV)" << endl;
  int size_ch = vec_2f3.size();
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_2f3[ch].size();
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int ap_Tz = state_T2P_Tz[ch][i].a;
      int bp_Tz = state_T2P_Tz[ch][i].b;

      int t2a_p = SP.state_Tz[ap_Tz].tz_2;
      int t2b_p = SP.state_Tz[bp_Tz].tz_2;
      int J2ab_p = state_T2P_Tz[ch][i].J_2;
      int T2ab_pz = state_T2P_Tz[ch][i].T_2z;
      int ap = std::floor(ap_Tz / 2.0);
      int bp = std::floor(bp_Tz / 2.0);

      for (int j = 0; j < size; j++) {
        int a_Tz = state_T2P_Tz[ch][j].a;
        int b_Tz = state_T2P_Tz[ch][j].b;
        int t2a = SP.state_Tz[a_Tz].tz_2;
        int t2b = SP.state_Tz[b_Tz].tz_2;
        int J2ab = state_T2P_Tz[ch][j].J_2;
        int T2ab_z = state_T2P_Tz[ch][j].T_2z;
        int a = std::floor(a_Tz / 2.0);
        int b = std::floor(b_Tz / 2.0);

        // if (T2ab_z == 0 && ap == bp && ap_Tz > 0)
        //  continue;
        // if (T2ab_z == 0 && a == b && a_Tz > 0)
        //  continue;
        file << ap_Tz + 1 << "\t" << bp_Tz + 1 << "\t" << a_Tz + 1 << "\t"
             << b_Tz + 1 << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3_C[ch][i][j] << endl;
        /*file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz <<
        "\t "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;
        /*
        file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
        1
             << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";
        file << vec_2f3[ch][i][j];
        file << "\t i = " << i << "\t j=" << j << "\t " << ap_Tz << "\t
        "
             << bp_Tz << "\t" << a_Tz << "\t" << b_Tz << "\t " << J2ab_p
             << "\t " << J2ab;
        file << endl;/**/
        // if (std::abs(vec_2f3[ch][i][j]) > 0.0001)
        // file << i << "\t" << j << "\t" << vec_2f3[ch][i][j] << endl;
      }
    }
  }
  file.close();
} /**/

void Norm_Ord::checkVsTok_pack_Tz(std::string filename1,
                                  std::string filename2) {
  std::ofstream file2;
  file2.open(filename2);

  std::ifstream file1;
  file1.open(filename1);
  char u1[256]; // unuse
  file1.getline(u1, 256, '\n');
  while (!file1.eof()) {
    int a, b, c, d, e, f;
    int a_x, b_x, c_x, d_x, e_x, f_x;
    int tz_a, tz_b, tz_c, tz_d, tz_e, tz_f;
    int J12, J12p, J2ab, J2de, J_2, T_2z;
    double J, val;
    file1 >> a >> b >> c >> J12 >> d >> e >> f >> J12p >> J >> tz_a >> tz_b >>
        tz_c >> tz_d >> tz_e >> tz_f >> val;

    J2ab = J12 * 2;
    J2de = J12p * 2;
    J_2 = std::floor(J * 2);
    T_2z = (tz_a + tz_b + tz_c) * 2 - 9;

    a_x = (a - 1) * 2;
    if (tz_a == 2)
      a_x += 1;

    b_x = (b - 1) * 2;
    if (tz_b == 2)
      b_x += 1;

    c_x = (c - 1) * 2;
    if (tz_c == 2)
      c_x += 1;

    d_x = (d - 1) * 2;
    if (tz_d == 2)
      d_x += 1;

    e_x = (e - 1) * 2;
    if (tz_e == 2)
      e_x += 1;

    f_x = (f - 1) * 2;
    if (tz_f == 2)
      f_x += 1;

    /*  cout << a_x << "\t" << b_x << "\t" << c_x << "\t" << J2ab <<
       "\t" << d_x
           << "\t" << d_x << "\t" << f_x << "\t" << J2de << "\t" << J_2
       <<
       "\t"
           << T_2z << "\t" << val << endl;/**/

    int x = TP.FindTP_Tz(J_2, T_2z, J2ab, a_x, b_x, c_x);
    int y = TP.FindTP_Tz(J_2, T_2z, J2de, d_x, e_x, f_x);

    int ch_Tz = TP.FindC_Tz(J_2, T_2z);
    if (this->vec_flag_Tz_S.val(ch_Tz, x, y) != 1) {
      cout << "Wrong happened @ Norm_Ord::checkVsTok_pack_Tz" << endl;
      cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
      exit(0);
    }
    // double val_z = this->vec_Tz[ch_Tz][x][y];
    double val_z = this->vec_Tz_S.val(ch_Tz, x, y);

    double diff = val - val_z;
    if (std::abs(diff) < 0.000001)
      diff = 0;

    file2 << a_x << "\t" << b_x << "\t" << c_x << "\t" << J2ab << "\t" << d_x
          << "\t" << d_x << "\t" << f_x << "\t" << J2de << "\t" << J_2 << "\t"
          << T_2z << "\t" << val << "\t "
          << "\t" << ch_Tz << "\t" << x << "\t" << y << "\t" << val_z << "\t "
          << diff << endl;
  }
  file1.close();
  file2.close();
}

void Norm_Ord::checkVsTok_norm2f3_Tz(std::string filename1,
                                     std::string filename2) {
  std::ofstream file2;
  file2.open(filename2);

  std::ifstream file1;
  file1.open(filename1);
  char u1[256]; // unuse
  file1.getline(u1, 256, '\n');
  file1.getline(u1, 256, '\n');
  file1.getline(u1, 256, '\n');
  while (!file1.eof()) {
    int a, b, d, e;
    int a_x, b_x, d_x, e_x;
    int J12, J2ab, tz, Tab_z2;
    double val_C, val_P, val;
    file1 >> a >> b >> d >> e >> J12 >> tz;
    file1 >> val_C;
    file1 >> val_P;
    J2ab = J12 * 2;
    Tab_z2 = 2 * tz;
    if (tz == 1) {
      // pp
      a_x = (a - 1) * 2 + 1;
      b_x = (b - 1) * 2 + 1;
      d_x = (d - 1) * 2 + 1;
      e_x = (e - 1) * 2 + 1;
    }
    if (tz == -1) {
      // nn
      a_x = (a - 1) * 2;
      b_x = (b - 1) * 2;
      d_x = (d - 1) * 2;
      e_x = (e - 1) * 2;
    }
    if (tz == 0) {
      // pn
      a_x = (a - 1) * 2 + 1;
      b_x = (b - 1) * 2;
      d_x = (d - 1) * 2 + 1;
      e_x = (e - 1) * 2;
    }

    /*cout << a << "\t" << b << "\t" << d << "\t" << e << "\t" << J12 <<
       "\t"
         << tz << "\t" << std::setw(6) << val << endl;/**/
    //<< "\t" << x << "\t" << y
    //<< std::setw(6) << val_z << "\t" << std::setw(6) << diff << endl;

    int ch_Tz = FindC_T2P_Tz(J2ab, Tab_z2);
    int x = Find_T2P_Tz(J2ab, Tab_z2, a_x, b_x);
    int y = Find_T2P_Tz(J2ab, Tab_z2, d_x, e_x);

    double val_z = vec_2f3[ch_Tz][x][y];
    double val_z_C = vec_2f3_C[ch_Tz][x][y];
    double val_z_P = vec_2f3_P[ch_Tz][x][y];

    double diff = val - val_z;
    if (std::abs(diff) < 0.0001)
      diff = 0.0;
    double diff_C = val_C - val_z_C;
    if (std::abs(diff_C) < 0.0001)
      diff_C = 0.0;
    double diff_P = val_P - val_z_P;
    if (std::abs(diff_P) < 0.0001)
      diff_P = 0.0;

    file2 << a << "\t " << b << "\t " << d << "\t " << e << J12 << "\t " << tz
          << "\t ";
    file2 << std::setw(6) << val_C << "\t" << std::setw(6) << val_z_C << "\t"
          << std::setw(6) << diff_C << "\t";
    file2 << std::setw(6) << val_P << "\t" << std::setw(6) << val_z_P << "\t"
          << std::setw(6) << diff_P << "\t";
    file2 << endl; /**/
    /*file2 << std::setw(6) << val_C << "\t " << x << "\t " << y <<
    "\t"; file2 << std::setw(6) << val_z_C << "\t" << std::setw(6) <<
    diff_C
          << endl; /**/
  }
  file1.close();
  file2.close();
}

void Norm_Ord::checkVsTok_norm2f3_Tz() {

  std::string filename1 = "data/Tok_data/norm_ord/spev_11BcDcE_check.dat";
  std::string filename2 = "data/Tok_data/norm_ord/normTz_TokVsYz.dat";

  std::ofstream file2;
  file2.open(filename2);
  file2 << "a\tb\td\te\tJ12\tTz\t val_C_Tok\t val_C_YZ\t diff\t\t "
           "val_P_Tok\t "
           "val_P_YZ\t diff"
        << endl;
  std::ifstream file1;
  file1.open(filename1);
  char u1[256]; // unuse
  file1.getline(u1, 256, '\n');
  file1.getline(u1, 256, '\n');
  file1.getline(u1, 256, '\n');
  while (!file1.eof()) {
    int a, b, d, e;
    int a_x, b_x, d_x, e_x;
    int J12, J2ab, tz, Tab_z2;
    double val_C, val_P, val;
    file1 >> a >> b >> d >> e >> J12 >> tz;
    file1 >> val_C;
    file1 >> val_P;
    J2ab = J12 * 2;
    Tab_z2 = 2 * tz;
    if (tz == 1) {
      // pp
      a_x = (a - 1) * 2 + 1;
      b_x = (b - 1) * 2 + 1;
      d_x = (d - 1) * 2 + 1;
      e_x = (e - 1) * 2 + 1;
    }
    if (tz == -1) {
      // nn
      a_x = (a - 1) * 2;
      b_x = (b - 1) * 2;
      d_x = (d - 1) * 2;
      e_x = (e - 1) * 2;
    }
    if (tz == 0) {
      // pn
      a_x = (a - 1) * 2 + 1;
      b_x = (b - 1) * 2;
      d_x = (d - 1) * 2 + 1;
      e_x = (e - 1) * 2;
    }

    /*cout << a << "\t" << b << "\t" << d << "\t" << e << "\t" << J12 <<
       "\t"
         << tz << "\t" << std::setw(6) << val_C << endl; /**/
    //<< "\t" << x << "\t" << y
    //<< std::setw(6) << val_z << "\t" << std::setw(6) << diff << endl;

    int ch_Tz = FindC_T2P_Tz(J2ab, Tab_z2);
    int x = Find_T2P_Tz(J2ab, Tab_z2, a_x, b_x);
    int y = Find_T2P_Tz(J2ab, Tab_z2, d_x, e_x);
    // cout<<"\t x = "<<x<<"\t y = "<<y<<"\t ch_Tz = "<<ch_Tz<<endl;
    double val_z = vec_2f3[ch_Tz][x][y];
    double val_z_C = vec_2f3_C[ch_Tz][x][y];
    double val_z_P = vec_2f3_P[ch_Tz][x][y];

    double diff = val - val_z;
    if (std::abs(diff) < 0.0001)
      diff = 0.0;
    double diff_C = val_C - val_z_C;
    if (std::abs(diff_C) < 0.0001)
      diff_C = 0.0;
    double diff_P = val_P - val_z_P;
    if (std::abs(diff_P) < 0.0001)
      diff_P = 0.0;

    file2 << a << "\t " << b << "\t " << d << "\t " << e << "\t" << J12 << "\t "
          << tz << "\t ";
    file2 << std::setw(6) << val_C << "\t  " << std::setw(6) << val_z_C << "\t "
          << std::setw(6) << diff_C << "\t\t ";
    file2 << std::setw(6) << val_P << "\t " << std::setw(6) << val_z_P << "\t "
          << std::setw(6) << diff_P << "\t ";
    file2 << endl; /**/
    /*file2 << std::setw(6) << val_C << "\t " << x << "\t " << y <<
    "\t"; file2 << std::setw(6) << val_z_C << "\t" << std::setw(6) <<
    diff_C
          << endl; /**/
  }
  file1.close();
  file2.close();
}

// =================== Normal Ord Sec ===================== //

void Norm_Ord_Sec::init(double hw, int E_cut, int sp_max, int sp_core,
                        int sp_P_beg, int sp_P_end, int sp_Q_beg,
                        int sp_Q_end) {
  this->hbar_omega = hw;
  cout << "!!! Norm_Ord_Sec Hbar_Omega = " << hw << endl;
  this->E_cut = E_cut;
  this->sp_max = sp_max;
  this->sp_core = sp_core;
  this->sp_P_beg = sp_P_beg;
  this->sp_P_end = sp_P_end;
  this->sp_Q_beg = sp_Q_beg;
  this->sp_Q_end = sp_Q_end;

  this->sp_max_Tz = sp_max * 2;
  this->sp_core_Tz = sp_core * 2 + 1;
  this->sp_P_beg_Tz = sp_P_beg * 2;
  this->sp_P_end_Tz = sp_P_end * 2 + 1;
  this->sp_Q_beg_Tz = sp_Q_beg * 2;
  this->sp_Q_end_Tz = sp_Q_end * 2 + 1;
  // init sp2rd
  int size_sp = SP.state_Tz.size();
  this->sp2rd.resize(size_sp);
  for (int i = 0; i < size_sp; i++) {
    int n = SP.state_Tz[i].n;
    int l = SP.state_Tz[i].l;
    int j_2 = SP.state_Tz[i].j_2;
    int tz_2 = SP.state_Tz[i].tz_2;
    int flag = -1;
    for (int j = 0; j < RD.sp_rd.state_Tz.size(); j++) {
      if (n != RD.sp_rd.state_Tz[j].n)
        continue;
      if (l != RD.sp_rd.state_Tz[j].l)
        continue;
      if (j_2 != RD.sp_rd.state_Tz[j].j_2)
        continue;
      if (tz_2 != (RD.sp_rd.state_Tz[j].tz_2 *
                   -1)) // !!! Morten: -1->proton Zhuo: -1->neutron
        continue;
      flag = j;
      break;
    }
    if (flag < 0) {
      cout << "Wrong @ Norm_Ord_Sec::init " << endl;
      cout << " this->sp2rd[i] = -1" << endl;
      exit(0);
    }
    this->sp2rd[i] = flag;
    // cout << "i : " << i << "\t sp2rd : " << flag << endl;
  }
}

void Norm_Ord_Sec::build_U() {
  // Totally U has three piece
  // 1: <h|u|p>  2: <h|u|a>  3: <a\U\p>
  int sp_size = SP.state_Tz.size();
  // cout << "sp_size : " << sp_size << endl;
  vec_U.resize(sp_size);
  vec_UT.resize(sp_size);
  vec_U_flag.resize(sp_size);
  vec_UT_flag.resize(sp_size);
  for (int i = 0; i < sp_size; i++) {
    vec_U[i].resize(sp_size);
    vec_UT[i].resize(sp_size);
    vec_U_flag[i].resize(sp_size);
    vec_UT_flag[i].resize(sp_size);
    for (int j = 0; j < sp_size; j++) {
      vec_U[i][j] = 0.0;
      vec_UT[i][j] = 0.0;
      vec_U_flag[i][j] = 0;
      vec_UT_flag[i][j] = 0;
    }
  }
  // 1: <h|u|p> 2: <h|u|a>
  for (int h = 0; h <= this->sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    for (int p = this->sp_P_beg_Tz; p <= this->sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      // cout << "\t -- h : " << h << "\t p : " << p << " --- " << endl;
      int E_dom_int = std::abs(e_h - e_p);
      if (E_dom_int > this->E_cut)
        continue;
      double val_U = this->cal1f2_Tz_C(h, p);
      vec_U[h][p] = val_U;
      vec_U_flag[h][p] = 1;
      double val_UT = this->cal1f3_Tz_C(h, p);
      vec_UT[h][p] = val_UT;
      vec_UT_flag[h][p] = 1;
      // cout << "1: <h|u|p> 2: <h|u|a> h : " << h << "\t p : " << p
      //      << "\t U = " << vec_U[h][p] << "\t UT = " << vec_UT[h][p] << endl;
    }
  }
  // 3: <a|u|p>
  for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
    int e_a = SP.state_Tz[a].E;
    for (int p = this->sp_Q_beg_Tz; p <= this->sp_Q_end_Tz; p++) {
      cout << "\t -- a : " << a << "\t p : " << p << " --- " << endl;
      int e_p = SP.state_Tz[p].E;
      int E_dom_int = std::abs(e_a - e_p);
      if (E_dom_int > this->E_cut)
        continue;
      double val_U = this->cal1f2_Tz_C(a, p);
      vec_U[a][p] = val_U;
      vec_U_flag[a][p] = 1;
      double val_UT = this->cal1f3_Tz_C(a, p);
      vec_UT[a][p] = val_UT;
      vec_UT_flag[a][p] = 1;
      // cout << "3: <a|u|p>  a : " << a << "\t p : " << p
      //      << "\t U = " << vec_U[a][p] << "\t UT = " << vec_UT[a][p] << endl;
    }
  }
}

void Norm_Ord_Sec::build_2V() {
  // Totally U has three piece
  // 1: <h1h2|V3N|ap>  2: <ah|V3N|ap>  3: <ah\V3N\p1p2>
  int size_ch = DP.state_Tz.size();
  // cout << "DP channel : " << size_ch << endl;
  vec_2V2N.resize(size_ch);
  vec_2V3N.resize(size_ch);
  vec_2V2N_flag.resize(size_ch);
  vec_2V3N_flag.resize(size_ch);
  for (int ch = 0; ch < size_ch; ch++) {
    int size = DP.state_Tz[ch].size();
    vec_2V2N[ch].resize(size);
    vec_2V3N[ch].resize(size);
    vec_2V2N_flag[ch].resize(size);
    vec_2V3N_flag[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2V2N[ch][i].resize(size);
      vec_2V3N[ch][i].resize(size);
      vec_2V2N_flag[ch][i].resize(size);
      vec_2V3N_flag[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2V2N[ch][i][j] = 0.0;
        vec_2V3N[ch][i][j] = 0.0;
        vec_2V2N_flag[ch][i][j] = 0;
        vec_2V3N_flag[ch][i][j] = 0;
      }
    }
  }
  // 1: <h1h2|V3N|ap>
  for (int h1 = 0; h1 <= this->sp_core_Tz; h1++) {
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    int j2_h1 = SP.state_Tz[h1].j_2;
    int l_h1 = SP.state_Tz[h1].l;
    int e_h1 = SP.state_Tz[h1].E;
    for (int h2 = 0; h2 <= this->sp_core_Tz; h2++) {
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      int j2_h2 = SP.state_Tz[h2].j_2;
      int l_h2 = SP.state_Tz[h2].l;
      int e_h2 = SP.state_Tz[h2].E;

      int T2z_x = t2z_h1 + t2z_h2;
      int Par_x = (l_h1 + l_h2) % 2;
      int J2_min = std::abs(j2_h1 - j2_h2);
      int J2_max = j2_h1 + j2_h2;

      for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
        int t2z_a = SP.state_Tz[a].tz_2;
        int j2_a = SP.state_Tz[a].j_2;
        int l_a = SP.state_Tz[a].l;
        int e_a = SP.state_Tz[a].E;
        for (int p = this->sp_P_beg_Tz; p <= this->sp_Q_end_Tz; p++) {
          int t2z_p = SP.state_Tz[p].tz_2;
          int j2_p = SP.state_Tz[p].j_2;
          int l_p = SP.state_Tz[p].l;
          int e_p = SP.state_Tz[p].E;

          if (T2z_x != (t2z_a + t2z_p))
            continue;
          int E_dom_int = std::abs(e_h1 + e_h2 - e_a - e_p);
          if (E_dom_int > this->E_cut)
            continue;
          // if (Par_x != (l_a + l_p) % 2)
          //   continue;

          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_a, j2_p) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, h1, h2);

            int index_y = DP.FindDP_Tz(J2, T2z_x, a, p);

            int h1_x = sp2rd[h1] + 1;
            int h2_x = sp2rd[h2] + 1;
            int a_x = sp2rd[a] + 1;
            int p_x = sp2rd[p] + 1;

            double mat = RD.get_tbme_unnormalized(h1_x, h2_x, a_x, p_x, J2);

            vec_2V2N[ch_Tz][index_x][index_y] = mat;
            vec_2V2N_flag[ch_Tz][index_x][index_y] = 1;

            double mat_3 = this->cal2f3_Tz_C(DP.state_Tz[ch_Tz][index_x],
                                             DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
          }
        }
      }
    }
  }

  // 2: <ah|V3N|ap>
  for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
    int t2z_a = SP.state_Tz[a].tz_2;
    int j2_a = SP.state_Tz[a].j_2;
    int l_a = SP.state_Tz[a].l;
    for (int h = 0; h <= this->sp_core_Tz; h++) {
      int t2z_h = SP.state_Tz[h].tz_2;
      int j2_h = SP.state_Tz[h].j_2;
      int l_h = SP.state_Tz[h].l;
      int e_h = SP.state_Tz[h].E;

      int T2z_x = t2z_a + t2z_h;
      int Par_x = (l_a + l_h) % 2;
      int J2_min = std::abs(j2_a - j2_h);
      int J2_max = j2_a + j2_h;

      for (int p = this->sp_P_beg_Tz; p <= this->sp_Q_end_Tz; p++) {
        int t2z_p = SP.state_Tz[p].tz_2;
        int j2_p = SP.state_Tz[p].j_2;
        int l_p = SP.state_Tz[p].l;
        int e_p = SP.state_Tz[p].E;
        if (t2z_h != t2z_p)
          continue;
        int E_dom_int = std::abs(e_h - e_p);
        if (E_dom_int > this->E_cut)
          continue;
        // if (Par_x != (l_a + l_p) % 2)
        //   continue;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          if (this->Tri(J2, j2_a, j2_p) < 0)
            continue;

          int ch_Tz = DP.FindC_Tz(J2, T2z_x);
          int index_x = DP.FindDP_Tz(J2, T2z_x, a, h);
          int index_y = DP.FindDP_Tz(J2, T2z_x, a, p);

          int a_x = sp2rd[a] + 1;
          int h_x = sp2rd[h] + 1;
          // int a_x = sp2rd[a] + 1;
          int p_x = sp2rd[p] + 1;

          double mat = RD.get_tbme_unnormalized(a_x, h_x, a_x, p_x, J2);

          vec_2V2N[ch_Tz][index_x][index_y] = mat;
          vec_2V2N_flag[ch_Tz][index_x][index_y] = 1;

          double mat_3 = this->cal2f3_Tz_C(DP.state_Tz[ch_Tz][index_x],
                                           DP.state_Tz[ch_Tz][index_y]);
          vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
          vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
        }
      }
    }
  }

  // 3: <ah|V3N|p1p2>
  for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
    int t2z_a = SP.state_Tz[a].tz_2;
    int j2_a = SP.state_Tz[a].j_2;
    int l_a = SP.state_Tz[a].l;
    int e_a = SP.state_Tz[a].E;
    for (int h = 0; h <= this->sp_core_Tz; h++) {
      int t2z_h = SP.state_Tz[h].tz_2;
      int j2_h = SP.state_Tz[h].j_2;
      int l_h = SP.state_Tz[h].l;
      int e_h = SP.state_Tz[h].E;

      int T2z_x = t2z_a + t2z_h;
      int Par_x = (l_a + l_h) % 2;
      int J2_min = std::abs(j2_a - j2_h);
      int J2_max = j2_a + j2_h;

      for (int p1 = this->sp_P_beg_Tz; p1 <= this->sp_Q_end_Tz; p1++) {
        int t2z_p1 = SP.state_Tz[p1].tz_2;
        int j2_p1 = SP.state_Tz[p1].j_2;
        int l_p1 = SP.state_Tz[p1].l;
        int e_p1 = SP.state_Tz[p1].E;
        for (int p2 = this->sp_P_beg_Tz; p2 <= this->sp_Q_end_Tz; p2++) {
          int t2z_p2 = SP.state_Tz[p2].tz_2;
          int j2_p2 = SP.state_Tz[p2].j_2;
          int l_p2 = SP.state_Tz[p2].l;
          int e_p2 = SP.state_Tz[p2].E;
          if (T2z_x != (t2z_p1 + t2z_p2))
            continue;
          if (Par_x != (l_p1 + l_p2) % 2)
            continue;
          int E_dom_int = std::abs(e_a + e_h - e_p1 - e_p2);
          if (E_dom_int > this->E_cut)
            continue;
          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_p1, j2_p2) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, a, h);
            int index_y = DP.FindDP_Tz(J2, T2z_x, p1, p2);

            int a_x = sp2rd[a] + 1;
            int h_x = sp2rd[h] + 1;
            int p1_x = sp2rd[p1] + 1;
            int p2_x = sp2rd[p2] + 1;

            double mat = RD.get_tbme_unnormalized(a_x, h_x, p1_x, p2_x, J2);

            vec_2V2N[ch_Tz][index_x][index_y] = mat;
            vec_2V2N_flag[ch_Tz][index_x][index_y] = 1;

            double mat_3 = this->cal2f3_Tz_C(DP.state_Tz[ch_Tz][index_x],
                                             DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
          }
        }
      }
    }
  } /**/
}
void Norm_Ord_Sec::build_2V_Luigi() {
  // Totally U has three piece
  //   1: <ai|V3N|jk>   2: <ib|V3N|jk>
  // (3): <jk|V3N|ai> (4): <jk|V3N|ib>
  int size_ch = DP.state_Tz.size();
  int E_cut = DP.E_cut_Luigi;
  int E_model_SP = DP.E_model_SP_Luigi;
  // cout << "DP channel : " << size_ch << endl;
  vec_2V2N.resize(size_ch);
  vec_2V2N_unnorm.resize(size_ch);

  vec_2V3N.resize(size_ch);
  vec_2V2N_flag.resize(size_ch);
  vec_2V3N_flag.resize(size_ch);
  for (int ch = 0; ch < size_ch; ch++) {
    int size = DP.state_Tz[ch].size();
    vec_2V2N[ch].resize(size);
    vec_2V2N_unnorm[ch].resize(size);
    vec_2V3N[ch].resize(size);
    vec_2V2N_flag[ch].resize(size);
    vec_2V3N_flag[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2V2N[ch][i].resize(size);
      vec_2V2N_unnorm[ch][i].resize(size);
      vec_2V3N[ch][i].resize(size);
      vec_2V2N_flag[ch][i].resize(size);
      vec_2V3N_flag[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2V2N[ch][i][j] = 0.0;
        vec_2V2N_unnorm[ch][i][j] = 0.0;
        vec_2V3N[ch][i][j] = 0.0;
        vec_2V2N_flag[ch][i][j] = 0;
        vec_2V3N_flag[ch][i][j] = 0;
      }
    }
  }
  // 1: <ai|V3N|jk>
  for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
    int t2z_a = SP.state_Tz[a].tz_2;
    int j2_a = SP.state_Tz[a].j_2;
    int l_a = SP.state_Tz[a].l;
    int e_a = SP.state_Tz[a].E;
    for (int i = 0; i <= this->sp_Q_end_Tz; i++) {
      int t2z_i = SP.state_Tz[i].tz_2;
      int j2_i = SP.state_Tz[i].j_2;
      int l_i = SP.state_Tz[i].l;
      int e_i = SP.state_Tz[i].E;

      int T2z_x = t2z_a + t2z_i;
      int Par_x = (l_a + l_i) % 2;
      int J2_min = std::abs(j2_a - j2_i);
      int J2_max = j2_a + j2_i;
      if (std::abs(2 * E_model_SP - e_a - e_i) > E_cut)
        continue;
      for (int j = 0; j <= this->sp_Q_end_Tz; j++) {
        int t2z_j = SP.state_Tz[j].tz_2;
        int j2_j = SP.state_Tz[j].j_2;
        int l_j = SP.state_Tz[j].l;
        int e_j = SP.state_Tz[j].E;
        for (int k = 0; k <= this->sp_Q_end_Tz; k++) {
          int t2z_k = SP.state_Tz[k].tz_2;
          int j2_k = SP.state_Tz[k].j_2;
          int l_k = SP.state_Tz[k].l;
          int e_k = SP.state_Tz[k].E;
          if (std::abs(2 * E_model_SP - e_k - e_j) > E_cut)
            continue;
          if (T2z_x != (t2z_j + t2z_k))
            continue;
          if (Par_x != (l_j + l_k) % 2)
            continue;

          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_j, j2_k) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, a, i);
            int index_y = DP.FindDP_Tz(J2, T2z_x, j, k);

            int a_x = sp2rd[a] + 1;
            int i_x = sp2rd[i] + 1;
            int j_x = sp2rd[j] + 1;
            int k_x = sp2rd[k] + 1;

            double mat = RD.get_tbme(a_x, i_x, j_x, k_x, J2);

            vec_2V2N[ch_Tz][index_x][index_y] = mat;
            vec_2V2N[ch_Tz][index_y][index_x] = mat;

            double mat_unnorm =
                RD.get_tbme_unnormalized(a_x, i_x, j_x, k_x, J2);

            vec_2V2N_unnorm[ch_Tz][index_x][index_y] = mat_unnorm;
            vec_2V2N_unnorm[ch_Tz][index_y][index_x] = mat_unnorm;
            vec_2V2N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V2N_flag[ch_Tz][index_y][index_x] = 1;

            double mat_3 = this->cal2f3_Tz_C(DP.state_Tz[ch_Tz][index_x],
                                             DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;
          }
        }
      }
    }
  }
  // 1: <ib|V3N|jk>

  for (int i = 0; i <= this->sp_Q_end_Tz; i++) {
    int t2z_i = SP.state_Tz[i].tz_2;
    int j2_i = SP.state_Tz[i].j_2;
    int l_i = SP.state_Tz[i].l;
    int e_i = SP.state_Tz[i].E;
    for (int b = this->sp_P_beg_Tz; b <= this->sp_P_end_Tz; b++) {
      int t2z_b = SP.state_Tz[b].tz_2;
      int j2_b = SP.state_Tz[b].j_2;
      int l_b = SP.state_Tz[b].l;
      int e_b = SP.state_Tz[b].E;

      if (std::abs(2 * E_model_SP - e_i - e_b) > E_cut)
        continue;
      int T2z_x = t2z_b + t2z_i;
      int Par_x = (l_b + l_i) % 2;
      int J2_min = std::abs(j2_b - j2_i);
      int J2_max = j2_b + j2_i;

      for (int j = 0; j <= this->sp_Q_end_Tz; j++) {
        int t2z_j = SP.state_Tz[j].tz_2;
        int j2_j = SP.state_Tz[j].j_2;
        int l_j = SP.state_Tz[j].l;
        int e_j = SP.state_Tz[j].E;
        for (int k = 0; k <= this->sp_Q_end_Tz; k++) {
          int t2z_k = SP.state_Tz[k].tz_2;
          int j2_k = SP.state_Tz[k].j_2;
          int l_k = SP.state_Tz[k].l;
          int e_k = SP.state_Tz[k].E;

          if (std::abs(2 * E_model_SP - e_k - e_j) > E_cut)
            continue;
          if (T2z_x != (t2z_j + t2z_k))
            continue;
          if (Par_x != (l_j + l_k) % 2)
            continue;

          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_j, j2_k) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, i, b);
            int index_y = DP.FindDP_Tz(J2, T2z_x, j, k);

            int i_x = sp2rd[i] + 1;
            int b_x = sp2rd[b] + 1;
            int j_x = sp2rd[j] + 1;
            int k_x = sp2rd[k] + 1;

            double mat = RD.get_tbme(i_x, b_x, j_x, k_x, J2);
            double mat_unnorm =
                RD.get_tbme_unnormalized(i_x, b_x, j_x, k_x, J2);

            vec_2V2N[ch_Tz][index_x][index_y] = mat;
            vec_2V2N[ch_Tz][index_y][index_x] = mat;
            vec_2V2N_unnorm[ch_Tz][index_x][index_y] = mat_unnorm;
            vec_2V2N_unnorm[ch_Tz][index_y][index_x] = mat_unnorm;
            vec_2V2N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V2N_flag[ch_Tz][index_y][index_x] = 1;

            double mat_3 = this->cal2f3_Tz_C(DP.state_Tz[ch_Tz][index_x],
                                             DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;
          }
        }
      }
    }
  }
}

double Norm_Ord_Sec::cal1f2_Tz_C(int a, int b) {
  double res = 0.0;
  int j2_a = SP.state_Tz[a].j_2;
  int j2_b = SP.state_Tz[b].j_2;
  int l_a = SP.state_Tz[a].l;
  int l_b = SP.state_Tz[b].l;
  int tz_a = SP.state_Tz[a].tz_2;
  int tz_b = SP.state_Tz[b].tz_2;
  if (j2_a != j2_b || l_a != l_b || tz_a != tz_b)
    return res;

  double val_t = 0.0;
  for (int h = 0; h <= this->sp_core_Tz; h++) {
    int l_h = SP.state_Tz[h].l;
    int Par = (l_h + l_a) % 2;
    int j2_h = SP.state_Tz[h].j_2;
    int J2_max = j2_a + j2_h;
    int J2_min = std::abs(j2_a - j2_h);
    for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
      int a_x = sp2rd[a] + 1;
      int b_x = sp2rd[b] + 1;
      int h_x = sp2rd[h] + 1;

      // cout << "J2 : " << J2 << "\t Par " << Par << "\t " << a_x << "\t " <<
      // h_x
      //      << "\t " << b_x << "\t " << h_x << "\t : " << endl;
      double mat = RD.get_tbme_unnormalized(a_x, h_x, b_x, h_x, J2);
      // cout << "\t\t\t " << mat << endl;
      val_t += (J2 + 1) * mat;
    }
  }
  res = 1.0 / (j2_a + 1) * val_t;
  return res;
}

double Norm_Ord_Sec::cal1f3_Tz_C(int a, int b) {
  double res = 0.0;

  int j2_a = SP.state_Tz[a].j_2;
  int j2_b = SP.state_Tz[b].j_2;
  int l_a = SP.state_Tz[a].l;
  int l_b = SP.state_Tz[b].l;
  int t2z_a = SP.state_Tz[a].tz_2;
  int t2z_b = SP.state_Tz[b].tz_2;
  if (j2_a != j2_b || l_a != l_b || t2z_a != t2z_b)
    return res;

  //
  // int j2_p = sp_t.j_2;
  // int t2z_p = sp_t.tz_2;
  // int p = sp_t.index;

  for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      int J2ab_min = std::abs(j2_h1 - j2_h2);
      int J2ab_max = j2_h1 + j2_h2;
      int T2z_x = t2z_a + t2z_h1 + t2z_h2;
      int T2z_y = t2z_b + t2z_h1 + t2z_h2;
      if (T2z_x != T2z_y)
        continue;
      int T2z = T2z_x;
      for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
        int J2_min = std::abs(J2ab - j2_a);
        int J2_max = J2ab + j2_a;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          int ch_Tz = TP.FindC_Tz(J2, T2z);
          int paul_cons_x = this->pauli_exc_Tz(h1, h2, a);
          int paul_cons_y = this->pauli_exc_Tz(h1, h2, b);

          if (paul_cons_x < 0 || paul_cons_y < 0)
            continue;
          int x = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, a);
          int y = TP.FindTP_Tz(J2, T2z, J2ab, h1, h2, b);
          // if (this->vec_3V_flag_Tz[ch_Tz][x][y] != 1) {
          if (this->vec_3V_flag_Tz_S.val(ch_Tz, x, y) != 1) {
            cout << "Wrong happened @ Norm_Ord_Sec::cal1f3_Tz" << endl;
            cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
            cout << "\t T2z : " << T2z << "\t J2 : " << J2
                 << "\t ch_Tz : " << ch_Tz << endl;
            cout << "\t ch_Tz : " << ch_Tz << "\t x : " << x << "\t y : " << y
                 << endl;
            exit(0);
          }
          // double mat_t = this->vec_3V_Tz[ch_Tz][x][y];
          double mat_t = this->vec_3V_Tz_S.val(ch_Tz, x, y);
          double C1 = (J2 + 1.0) / (2.0 * (j2_a + 1.0));
          double mat = C1 * mat_t;

          // cout << "\t J2 = " << J2 << "\t T2z = " << T2z
          //     << "\t J2ab = " << J2ab<<"\t j2_p = "<<j2_p<<"\t";
          //   << endl;
          // cout << "h1 = " << h1 << "\t h2 = " << h2 << "\t p = " << p
          // cout<<"C1 = "<<C1 << "\t mat_t = " << mat_t <<"\t mat = "<<
          // mat << endl; /**/
          res += mat;
        }
      }
    }
  }
  // res = res;
  return res;
}

double Norm_Ord_Sec::cal2f3_Tz_C(State_DP &T2P_i, State_DP &T2P_f) {
  double res = 0.0;

  int T2ab_z_i = T2P_i.T_2z;
  int T2ab_z_f = T2P_f.T_2z;
  int J2ab_i = T2P_i.J_2;
  int J2ab_f = T2P_f.J_2;
  if (T2ab_z_i != T2ab_z_f)
    return res;
  if (J2ab_i != J2ab_f)
    return res;
  int T2ab_z = T2ab_z_i;
  int J2ab = J2ab_i;
  int a = T2P_i.a;
  int b = T2P_i.b;
  int d = T2P_f.a;
  int e = T2P_f.b;

  int h_core = this->sp_core_Tz;

  // cout << "--- -- a = " << a << "\t b = " << b << "\t\t d = " << d
  //     << "\t e = " << e << endl;

  // double result = 0.0;
  for (int h = 0; h <= h_core; h++) {
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2max = j2_h + J2ab;
    int J2min = std::abs(j2_h - J2ab);
    int T2z = T2ab_z + t2z_h;
    for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) {
      int ch_Tz = TP.FindC_Tz(J2, T2z);
      /*cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " <<
      J2ab
           << endl;
      cout << "a = " << a << "\t b = " << b << "\t h = " << h <<
      endl;/**/
      if (SP.state_Tz[a].l == 0 && a == b && b == h)
        continue;
      if (SP.state_Tz[d].l == 0 && d == e && e == h)
        continue; /**/
      int x = TP.FindTP_Tz(J2, T2z, J2ab, a, b, h);
      // cout << "d = " << d << "\t e = " << e << "\t h = " << h <<
      // endl;/**/
      int y = TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
      if (x < 0) {
        cout << "Wrong at cal2f3_Tz x< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
        exit(0);
      }
      if (y < 0) {
        cout << "Wrong at cal2f3_Tz y< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "d = " << d << "\t e = " << e << "\t h = " << h << endl;
        exit(0);
      }
      // cout << "ch_Tz = " << ch_Tz << "\t x = " << x << "\t y = " << y
      // << endl;
      // if (this->vec_3V_flag_Tz[ch_Tz][x][y] != 1) {
      if (this->vec_3V_flag_Tz_S.val(ch_Tz, x, y) != 1) {
        cout << "Wrong happened @ Norm_Ord::cal_val2f3_Tz 1" << endl;
        cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
        exit(0);
      }
      // double mat_t = this->vec_3V_Tz[ch_Tz][x][y];
      double mat_t = this->vec_3V_Tz_S.val(ch_Tz, x, y);
      // cout << "mat_t = " << mat_t << endl;
      double mat = (J2 + 1.0) / (J2ab + 1.0) * mat_t;
      res += mat;
    }
  }
  return res;
}

void Norm_Ord_Sec::build_sec() {

  cout << " this->hbar_omega = " << this->hbar_omega << endl;
  int sp_size = SP.state_Tz.size();
  vec_A1.resize(sp_size);
  vec_A2.resize(sp_size);
  vec_B1.resize(sp_size);
  vec_B2.resize(sp_size);
  vec_C1.resize(sp_size);
  vec_C2.resize(sp_size);
  vec_D1.resize(sp_size);
  vec_D2.resize(sp_size);
  vec_E1.resize(sp_size);
  vec_E2.resize(sp_size);
  vec_E3.resize(sp_size);
  vec_F1.resize(sp_size);
  vec_F2.resize(sp_size);
  vec.resize(sp_size);
  vec_sub.resize(13);
  for (int i = 0; i < 13; i++) {
    vec_sub[i].resize(sp_size);
    for (int j = 0; j < sp_size; j++) {
      vec_sub[i][j] = 0.0;
    }
  }
  for (int i = 0; i < sp_size; i++) {
    vec_A1[i] = 0.0;
    vec_A2[i] = 0.0;
    vec_B1[i] = 0.0;
    vec_B2[i] = 0.0;
    vec_C1[i] = 0.0;
    vec_C2[i] = 0.0;
    vec_D1[i] = 0.0;
    vec_D2[i] = 0.0;
    vec_E1[i] = 0.0;
    vec_E2[i] = 0.0;
    vec_E3[i] = 0.0;
    vec_F1[i] = 0.0;
    vec_F2[i] = 0.0;
    vec[i] = 0.0;
  }
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    vec_A1[a] = cal_A1(a);
    vec_A2[a] = cal_A2(a);
    vec_B1[a] = cal_B1(a);
    vec_B2[a] = cal_B2(a);
    vec_C1[a] = cal_C1(a);
    vec_C2[a] = cal_C2(a);
    vec_D1[a] = cal_D1(a);
    vec_D2[a] = cal_D2(a);
    vec_E1[a] = cal_E1(a);
    vec_E2[a] = cal_E2(a);
    vec_E3[a] = cal_E3(a);
    vec_F1[a] = 0.0; // cal_F1(a);
    vec_F2[a] = 0.0; // cal_F2(a);
    // vec[a] = vec_A1[a] + vec_A2[a] + vec_B1[a] + vec_B2[a] + vec_C1[a] +
    //          vec_C2[a] + vec_D1[a] + vec_D2[a] + vec_E1[a] + vec_E2[a] +
    //          vec_E3[a] + vec_F1[a] + vec_F2[a];
    vec[a] = vec_A1[a] * 2 + vec_A2[a] * 2 + vec_B1[a] + vec_B2[a] +
             vec_C1[a] * 2 + vec_C2[a] * 2 + vec_D1[a] + vec_D2[a] +
             vec_E1[a] * 2 + vec_E2[a] * 2 + vec_E3[a] * 2 + vec_F1[a] +
             vec_F2[a];
    vec_sub[0][a] = vec_A1[a];
    vec_sub[1][a] = vec_A2[a];
    vec_sub[2][a] = vec_B1[a];
    vec_sub[3][a] = vec_B2[a];
    vec_sub[4][a] = vec_C1[a];
    vec_sub[5][a] = vec_C2[a];
    vec_sub[6][a] = vec_D1[a];
    vec_sub[7][a] = vec_D2[a];
    vec_sub[8][a] = vec_E1[a];
    vec_sub[9][a] = vec_E2[a];
    vec_sub[10][a] = vec_E3[a];
    vec_sub[11][a] = vec_F1[a];
    vec_sub[12][a] = vec_F2[a];
    // cout << "a : " << a << "\t " << vec_A1[a] << "\t " << vec_A2[a] << "\t
    // "
    //      << vec_B1[a] << "\t " << vec_B2[a] << "\t " << vec_C1[a] << "\t "
    //      << vec_C2[a] << "\t " << vec_D1[a] << "\t " << vec_D2[a] << "\t "
    //      << vec_E1[a] << "\t " << vec_E2[a] << "\t " << vec_E3[a] << "\t "
    //      << vec_F1[a] << "\t " << vec_F2[a] << endl;
  }
  cout << "a\t A1*2\t A2*2\t B1\t B2\t C1*2\t C2*2\t D1\t D2\t E1*2\t E2*2\t "
          "E3*2\t F1\t "
          "F2\t  "
       << endl;
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    cout << a << "\t " << std::setw(8) << std::setprecision(5) << vec_A1[a] * 2
         << "\t " << std::setw(8) << std::setprecision(5) << vec_A2[a] * 2
         << "\t " << std::setw(8) << std::setprecision(5) << vec_B1[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_B2[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_C1[a] * 2 << "\t "
         << std::setw(8) << std::setprecision(5) << vec_C2[a] * 2 << "\t "
         << std::setw(8) << std::setprecision(5) << vec_D1[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_D2[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_E1[a] * 2 << "\t "
         << std::setw(8) << std::setprecision(5) << vec_E2[a] * 2 << "\t "
         << std::setw(8) << std::setprecision(5) << vec_E3[a] * 2 << "\t "
         << std::setw(8) << std::setprecision(5) << vec_F1[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_F2[a] << endl;
  }
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    cout << " a : " << a << "\t vec_D1[a] : " << vec_D1[a]
         << "\t vec_E1[a] : " << vec_E1[a] << endl;
  }
}

double Norm_Ord_Sec::cal_A1(int a) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    double E_dom = this->hbar_omega * (e_a - e_p);
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (E_dom == 0.0)
      continue;
    if (j2_a != j2_p || t2z_a != t2z_p)
      continue;
    if (vec_U_flag[a][p] != 1 || vec_UT_flag[a][p] != 1) {
      cout << "Wrong @ Norm_Ord_Sec::cal_A1" << endl;
      cout << "vec_U_flag[a][p] : " << vec_U_flag[a][p];
      cout << "\t vec_UT_flag[a][p] : " << vec_UT_flag[a][p] << endl;
      cout << "a : " << a << "\t p : " << p << endl;
      exit(0);
    }
    double res_t = vec_U[a][p] * vec_UT[a][p];
    res_t = res_t / E_dom;
    res += res_t;
  }
  double result = res * phase;
  return result;
}

double Norm_Ord_Sec::cal_A2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    double E_dom = this->hbar_omega * (e_h - e_a);
    if (E_dom == 0.0)
      continue;
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (j2_a != j2_h || t2z_a != t2z_h)
      continue;
    if (vec_U_flag[h][a] != 1 || vec_UT_flag[h][a] != 1) {
      cout << "Wrong @ Norm_Ord_Sec::cal_A2" << endl;
      cout << "vec_U_flag[h][a] : " << vec_U_flag[h][a];
      cout << "\t vec_UT_flag[h][a] : " << vec_UT_flag[h][a] << endl;
      cout << "\t h : " << h << "\t a : " << a << endl;
      exit(0);
    }
    double res_t = vec_U[h][a] * vec_UT[h][a];
    res_t = res_t / E_dom;
    res += res_t;
  }
  double result = res * phase;
  return result;
}

double Norm_Ord_Sec::cal_B1(int a) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    double E_dom = this->hbar_omega * (e_a - e_p);
    if (E_dom == 0.0)
      continue;
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (j2_a != j2_p || t2z_a != t2z_p)
      continue;
    if (vec_UT_flag[a][p] != 1) {
      cout << "Wrong @ Norm_Ord_Sec::cal_B1" << endl;
      cout << "vec_UT_flag[a][p] : " << vec_UT_flag[a][p] << endl;
      cout << "a : " << a << "\t p : " << p << endl;
      exit(0);
    }
    double res_t = vec_UT[a][p] * vec_UT[a][p];
    res_t = res_t / E_dom;
    res += res_t;
  }
  double result = res * phase;
  return result;
}

double Norm_Ord_Sec::cal_B2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    double E_dom = this->hbar_omega * (e_h - e_a);
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (E_dom == 0.0)
      continue;
    if (j2_a != j2_h || t2z_a != t2z_h)
      continue;
    if (vec_U_flag[h][a] != 1 || vec_UT_flag[h][a] != 1) {
      cout << "Wrong @ Norm_Ord_Sec::cal_B2" << endl;
      cout << "\t vec_UT_flag[h][a] : " << vec_UT_flag[h][a] << endl;
      cout << "\t h : " << h << "\t a : " << a << endl;
      exit(0);
    }
    double res_t = vec_UT[h][a] * vec_UT[h][a];
    res_t = res_t / E_dom;
    res += res_t;
  }
  double result = res * phase;
  return result;
}

double Norm_Ord_Sec::cal_C1(int a) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2_min = std::abs(j2_a - j2_h);
    int J2_max = j2_a + j2_h;
    int T2z = t2z_a + t2z_h;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[h].l;
    int Par_x = par_t % 2;

    for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
      int e_p1 = SP.state_Tz[p1].E;
      int j2_p1 = SP.state_Tz[p1].j_2;
      int t2z_p1 = SP.state_Tz[p1].tz_2;
      for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
        int e_p2 = SP.state_Tz[p2].E;
        int j2_p2 = SP.state_Tz[p2].j_2;
        int t2z_p2 = SP.state_Tz[p2].tz_2;
        if (T2z != (t2z_p1 + t2z_p2))
          continue;
        par_t = SP.state_Tz[p1].l + SP.state_Tz[p2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          // cout << "-- C1 --  J2 : " << J2 << "\t j2_p1 : " << j2_p1
          //      << "\t j2_p2 : " << j2_p2 << "\t a : " << a << "\t h : " <<
          //      h
          //      << "\t p1 : " << p1 << "\t p2 : " << p2 << endl;
          if (Tri(J2, j2_p1, j2_p2) < 0)
            continue;
          double E_dom = this->hbar_omega * (e_a + e_h - e_p1 - e_p2);
          if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
            continue;
          int ch_DP = DP.FindC_Tz(J2, T2z);
          // cout << "\t \t ch_DP = " << ch_DP << endl;

          int index_x = DP.FindDP_Tz(J2, T2z, a, h);
          int index_y = DP.FindDP_Tz(J2, T2z, p1, p2);
          if (vec_2V2N_flag[ch_DP][index_x][index_y] != 1 ||
              vec_2V3N_flag[ch_DP][index_x][index_y] != 1) {
            cout << "Wrong @ Norm_Ord_Sec::cal_C1" << endl;
            cout << "vec_2V2N_flag[ch_DP][x][y] : "
                 << vec_2V2N_flag[ch_DP][index_x][index_y];
            cout << "\t vec_2V3N_flag[ch_DP][x][y] : "
                 << vec_2V3N_flag[ch_DP][index_x][index_y] << endl;
            cout << "\t ch_DP : " << ch_DP << "\t index_x : " << index_x
                 << "\t index_y : " << index_y << endl;
            exit(0);
          }
          double res_t = vec_2V2N[ch_DP][index_x][index_y] *
                         vec_2V3N[ch_DP][index_x][index_y];

          res += (J2 + 1) * res_t / E_dom;
          // cout << "\t C1 += " << res_t << "\t";
          // cout << "a : " << a << "\t h : " << h << "\t p1 : " << p1
          //      << "\t p2 : " << p2 << endl;
        }
      }
    }
  }
  double result = res * phase / (2.0 * (j2_a + 1.0));
  return result;
}

double Norm_Ord_Sec::cal_C2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  double res = 0.0;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    int J2_min = std::abs(j2_a - j2_p);
    int J2_max = j2_a + j2_p;
    int T2z = t2z_a + t2z_p;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[p].l;
    int Par_x = par_t % 2;

    for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
      int e_h1 = SP.state_Tz[h1].E;
      int j2_h1 = SP.state_Tz[h1].j_2;
      int t2z_h1 = SP.state_Tz[h1].tz_2;
      for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
        int e_h2 = SP.state_Tz[h2].E;
        int j2_h2 = SP.state_Tz[h2].j_2;
        int t2z_h2 = SP.state_Tz[h2].tz_2;
        double E_dom = this->hbar_omega * (e_h1 + e_h2 - e_a - e_p);
        if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
          continue;
        if (T2z != (t2z_h1 + t2z_h2))
          continue;
        par_t = SP.state_Tz[h1].l + SP.state_Tz[h2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          // cout << "-- C2 --  J2 : " << J2 << "\t j2_h1 : " << j2_h1
          //      << "\t j2_h1 : " << j2_h2 << "\t a : " << a << "\t p : " <<
          //      p
          //      << "\t h1 : " << h1 << "\t h2 : " << h2 << endl;
          if (Tri(J2, j2_h1, j2_h2) < 0)
            continue;
          int ch_DP = DP.FindC_Tz(J2, T2z);
          // cout << "\t \t ch_DP = " << ch_DP << endl;
          int index_x = DP.FindDP_Tz(J2, T2z, a, p);

          int index_y = DP.FindDP_Tz(J2, T2z, h1, h2);
          if (vec_2V2N_flag[ch_DP][index_y][index_x] != 1 ||
              vec_2V3N_flag[ch_DP][index_y][index_x] != 1) {
            cout << "Wrong @ Norm_Ord_Sec::cal_C2" << endl;
            cout << "vec_2V2N_flag[ch_DP][y][x] : "
                 << vec_2V2N_flag[ch_DP][index_y][index_x];
            cout << "\t vec_2V3N_flag[ch_DP][y][x] : "
                 << vec_2V3N_flag[ch_DP][index_y][index_x] << endl;
            cout << "\t ch_DP : " << ch_DP << "\t index_y : " << index_y
                 << "\t index_x : " << index_x << endl;
            cout << "a : " << a << "\t p : " << p << "\t h1 : " << h1
                 << "\t h2 : " << h2 << endl;
            exit(0);
          }
          double res_t = vec_2V2N[ch_DP][index_y][index_x] *
                         vec_2V3N[ch_DP][index_y][index_x];

          res += (J2 + 1) * res_t / E_dom;
          // cout << "\t C2 += " << res_t << "\t";
          // cout << "a : " << a << "\t p : " << p << "\t h1 : " << h1
          //      << "\t h2 : " << h2 << endl;
        }
      }
    }
  }
  double result = res * phase / (2.0 * (j2_a + 1.0));
  return result;
}

double Norm_Ord_Sec::cal_D1(int a) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  // cout << "cal_D1 :  a = " << a << endl;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2_min = std::abs(j2_a - j2_h);
    int J2_max = j2_a + j2_h;
    int T2z = t2z_a + t2z_h;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[h].l;
    int Par_x = par_t % 2;

    for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
      int e_p1 = SP.state_Tz[p1].E;
      int j2_p1 = SP.state_Tz[p1].j_2;
      int t2z_p1 = SP.state_Tz[p1].tz_2;

      for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
        int e_p2 = SP.state_Tz[p2].E;
        int j2_p2 = SP.state_Tz[p2].j_2;
        int t2z_p2 = SP.state_Tz[p2].tz_2;
        par_t = SP.state_Tz[p1].l + SP.state_Tz[p2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        if (T2z != (t2z_p1 + t2z_p2))
          continue;
        double E_dom = this->hbar_omega * (e_a + e_h - e_p1 - e_p2);
        if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
          continue;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          if (Tri(J2, j2_p1, j2_p2) < 0)
            continue;
          int ch_DP = DP.FindC_Tz(J2, T2z);
          int index_x = DP.FindDP_Tz(J2, T2z, a, h);
          int index_y = DP.FindDP_Tz(J2, T2z, p1, p2);
          if (vec_2V3N_flag[ch_DP][index_x][index_y] != 1) {
            cout << "Wrong @ Norm_Ord_Sec::cal_D1" << endl;

            cout << "\t vec_2V3N_flag[ch_DP][x][y] : "
                 << vec_2V3N_flag[ch_DP][index_x][index_y] << endl;
            cout << "\t ch_DP : " << ch_DP << "\t index_x : " << index_x
                 << "\t index_y : " << index_y << endl;
            exit(0);
          }
          double res_t = vec_2V3N[ch_DP][index_x][index_y] *
                         vec_2V3N[ch_DP][index_x][index_y];

          double res_add = (J2 + 1) * res_t / E_dom;
          res += res_add;
          cout << "cal_D1 : J2 : " << J2 << "\t T2z : " << T2z << "\t a : " << a
               << "\t h : " << h << "\t p1 : " << p1 << "\t p2 : " << p2
               << endl;
          cout << "\t \t res : " << res << "\t E_dom : " << E_dom
               << "\t res_add : " << res_add << endl;
        }
      }
    }
  }
  double result = res * phase / (2.0 * (j2_a + 1.0));
  cout << "----- cal_D1 result = " << result << endl;
  return result;
}

double Norm_Ord_Sec::cal_D2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  double res = 0.0;
  // cout << "cal_D2 :  a = " << a;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    int J2_min = std::abs(j2_a - j2_p);
    int J2_max = j2_a + j2_p;
    int T2z = t2z_a + t2z_p;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[p].l;
    int Par_x = par_t % 2;
    // cout << "\t \t p = " << p << endl;

    for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
      int e_h1 = SP.state_Tz[h1].E;
      int j2_h1 = SP.state_Tz[h1].j_2;
      int t2z_h1 = SP.state_Tz[h1].tz_2;
      // cout << "\t h1 = " << h1;
      for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
        int e_h2 = SP.state_Tz[h2].E;
        int j2_h2 = SP.state_Tz[h2].j_2;
        int t2z_h2 = SP.state_Tz[h2].tz_2;
        if (T2z != (t2z_h1 + t2z_h2))
          continue;
        par_t = SP.state_Tz[h1].l + SP.state_Tz[h2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        double E_dom = this->hbar_omega * (e_h1 + e_h2 - e_a - e_p);
        if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
          continue;
        // cout << "\t h2 = " << h2;
        // cout << "\t E_dom : " << E_dom << endl;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          if (Tri(J2, j2_h1, j2_h2) < 0)
            continue;
          int ch_DP = DP.FindC_Tz(J2, T2z);
          int index_x = DP.FindDP_Tz(J2, T2z, a, p);
          int index_y = DP.FindDP_Tz(J2, T2z, h1, h2);
          if (vec_2V3N_flag[ch_DP][index_y][index_x] != 1) {
            cout << "Wrong @ Norm_Ord_Sec::cal_D2" << endl;
            cout << "\t vec_2V3N_flag[ch_DP][y][x] : "
                 << vec_2V3N_flag[ch_DP][index_y][index_x] << endl;
            cout << "\t ch_DP : " << ch_DP << "\t index_y : " << index_y
                 << "\t index_x : " << index_x << endl;
            exit(0);
          }
          double res_t = vec_2V3N[ch_DP][index_y][index_x] *
                         vec_2V3N[ch_DP][index_y][index_x];
          res_t = (J2 + 1) * res_t / E_dom;
          res += res_t;
          // cout << "\t \t ---- J2 = " << J2 << "\t res_t = " << res_t <<
          // endl;
        }
      }
    }
  }
  double result = res * phase / (2.0 * (j2_a + 1.0));
  // cout << "==== D2 result : " << result << endl;
  return result;
}

double Norm_Ord_Sec::cal_E1(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_ah = t2z_a + t2z_h;
    int J2_max = j2_a + j2_h;
    int J2_min = std::abs(j2_a - j2_h);
    int Par_h = SP.state_Tz[h].l % 2;
    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      int Par_p = SP.state_Tz[p].l % 2;
      double E_dom = this->hbar_omega * (e_h - e_p);
      if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
        continue;
      int T2z_ap = t2z_a + t2z_p;
      if (j2_h != j2_p || t2z_h != t2z_p || T2z_ap != T2z_ah)
        continue;
      if (Par_h != Par_p)
        continue;
      if (vec_U_flag[h][p] != 1) {
        cout << "Wrong @ Norm_Ord_Sec::cal_E1 : 1" << endl;
        cout << " vec_U_flag[h][p] : " << vec_U_flag[h][p] << endl;
        cout << "\t h : " << h << "\t p : " << p << endl;
        exit(0);
      }
      double val_U = vec_U[h][p];

      // double val_t = 0.0;
      for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
        if (Tri(J2, j2_a, j2_p) < 0)
          continue;
        int ch_DP = DP.FindC_Tz(J2, T2z_ap);
        int index_x = DP.FindDP_Tz(J2, T2z_ap, a, p);
        int index_y = DP.FindDP_Tz(J2, T2z_ah, a, h);
        if (vec_2V3N_flag[ch_DP][index_y][index_x] != 1) {
          cout << "Wrong @ Norm_Ord_Sec::cal_E1 : 2" << endl;
          cout << "\t vec_2V3N_flag[ch_DP][y][x] : "
               << vec_2V3N_flag[ch_DP][index_y][index_x] << endl;
          cout << "\t ch_DP : " << ch_DP << "\t index_y : " << index_y
               << "\t index_x : " << index_x << endl;
          exit(0);
        }
        res += (J2 + 1) * val_U * vec_2V3N[ch_DP][index_y][index_x] / E_dom;
      }
    }
  }
  double result = res * phase / (j2_a + 1.0);
  return result;
}

double Norm_Ord_Sec::cal_E2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int Par_h = SP.state_Tz[h].l % 2;
    int T2z_ah = t2z_a + t2z_h;
    int J2_max = j2_a + j2_h;
    int J2_min = std::abs(j2_a - j2_h);
    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      int Par_p = SP.state_Tz[p].l % 2;
      int T2z_ap = t2z_a + t2z_p;
      if (j2_h != j2_p || t2z_h != t2z_p || T2z_ap != T2z_ah)
        continue;
      if (Par_h != Par_p)
        continue;
      double E_dom = this->hbar_omega * (e_h - e_p);
      if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
        continue;
      if (vec_UT_flag[h][p] != 1) {
        cout << "Wrong @ Norm_Ord_Sec::cal_E2 : 1" << endl;
        cout << " vec_UT_flag[h][p] : " << vec_UT_flag[h][p] << endl;
        cout << "\t h : " << h << "\t p : " << p << endl;
        exit(0);
      }
      double val_UT = vec_UT[h][p];

      // double val_t = 0.0;
      for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
        if (Tri(J2, j2_a, j2_p) < 0)
          continue;
        int ch_DP = DP.FindC_Tz(J2, T2z_ap);
        int index_x = DP.FindDP_Tz(J2, T2z_ap, a, p);
        int index_y = DP.FindDP_Tz(J2, T2z_ah, a, h);
        if (vec_2V2N_flag[ch_DP][index_y][index_x] != 1) {
          cout << "Wrong @ Norm_Ord_Sec::cal_E2 : 2" << endl;
          cout << "\t vec_2V2N_flag[ch_DP][y][x] : "
               << vec_2V2N_flag[ch_DP][index_y][index_x] << endl;
          cout << "\t ch_DP : " << ch_DP << "\t index_y : " << index_y
               << "\t index_x : " << index_x << endl;
          exit(0);
        }
        res += (J2 + 1) * val_UT * vec_2V2N[ch_DP][index_y][index_x] / E_dom;
      }
    }
  }
  double result = res * phase / (j2_a + 1.0);
  return result;
}

double Norm_Ord_Sec::cal_E3(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int Par_h = SP.state_Tz[h].l % 2;
    int T2z_ah = t2z_a + t2z_h;
    int J2_max = j2_a + j2_h;
    int J2_min = std::abs(j2_a - j2_h);
    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      int Par_p = SP.state_Tz[p].l % 2;
      int T2z_ap = t2z_a + t2z_p;
      if (j2_h != j2_p || t2z_h != t2z_p || T2z_ap != T2z_ah)
        continue;
      if (Par_h != Par_p)
        continue;
      double E_dom = this->hbar_omega * (e_h - e_p);
      if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
        continue;
      if (vec_UT_flag[h][p] != 1) {
        cout << "Wrong @ Norm_Ord_Sec::cal_E3 : 1" << endl;
        cout << " vec_UT_flag[h][p] : " << vec_UT_flag[h][p] << endl;
        cout << "\t h : " << h << "\t p : " << p << endl;
        exit(0);
      }
      double val_UT = vec_UT[h][p];

      // /double val_t = 0.0;
      for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
        if (Tri(J2, j2_a, j2_p) < 0)
          continue;
        int ch_DP = DP.FindC_Tz(J2, T2z_ap);
        int index_x = DP.FindDP_Tz(J2, T2z_ap, a, p);
        int index_y = DP.FindDP_Tz(J2, T2z_ah, a, h);
        if (vec_2V3N_flag[ch_DP][index_y][index_x] != 1) {
          cout << "Wrong @ Norm_Ord_Sec::cal_E3 : 2" << endl;
          cout << "\t vec_2V3N_flag[ch_DP][y][x] : "
               << vec_2V3N_flag[ch_DP][index_y][index_x] << endl;
          cout << "\t ch_DP : " << ch_DP << "\t index_y : " << index_y
               << "\t index_x : " << index_x << endl;
          exit(0);
        }
        res += (J2 + 1) * val_UT * vec_2V3N[ch_DP][index_y][index_x] / E_dom;
      }
    }
  }
  double result = res * phase / (j2_a + 1.0);
  return result;
}

double Norm_Ord_Sec::cal_F1(int a) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;

  for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
    int e_h1 = SP.state_Tz[h1].E;
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    int J2_ah1_min = std::abs(j2_a - j2_h1);
    int J2_ah1_max = j2_a + j2_h1;
    for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
      int e_h2 = SP.state_Tz[h2].E;
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      if (pauli_exc_Tz(a, h1, h2) < 0)
        continue;
      int T2z = t2z_a + t2z_h1 + t2z_h2;

      for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
        int e_p = SP.state_Tz[p].E;
        int j2_p = SP.state_Tz[p].j_2;
        int t2z_p = SP.state_Tz[p].tz_2;
        for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
          int e_p1 = SP.state_Tz[p1].E;
          int j2_p1 = SP.state_Tz[p1].j_2;
          int t2z_p1 = SP.state_Tz[p1].tz_2;
          int J2_pp1_min = std::abs(j2_p - j2_p1);
          int J2_pp1_max = j2_p + j2_p1;
          for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
            int e_p2 = SP.state_Tz[p2].E;
            int j2_p2 = SP.state_Tz[p2].j_2;
            int t2z_p2 = SP.state_Tz[p2].tz_2;
            int T2z_y = t2z_p + t2z_p1 + t2z_p2;
            if (pauli_exc_Tz(p, p1, p2) < 0)
              continue;
            if (T2z_y != T2z)
              continue;

            double E_dom =
                this->hbar_omega * (e_a + e_h1 + e_h2 - e_p - e_p1 - e_p2);
            if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
              continue;
            for (int J2_ah1 = J2_ah1_min; J2_ah1 <= J2_ah1_max;
                 J2_ah1 = J2_ah1 + 2) {
              int J2_min = std::abs(J2_ah1 - j2_h2);
              int J2_max = J2_ah1 + j2_h2;
              for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
                for (int J2_pp1 = J2_pp1_min; J2_pp1 <= J2_pp1_max;
                     J2_pp1 = J2_pp1 + 2) {
                  if (Tri(J2, J2_pp1, j2_p2) < 0)
                    continue;
                  int ch_TP = TP.FindC_Tz(J2, T2z);
                  int index_x = TP.FindTP_Tz(J2, T2z, J2_ah1, a, h1, h2);
                  int index_y = TP.FindTP_Tz(J2, T2z, J2_pp1, p, p1, p2);
                  // if (vec_3V_flag_Tz[ch_TP][index_x][index_y] != 1) {
                  if (vec_3V_flag_Tz_S.val(ch_TP, index_x, index_y) != 1) {
                    cout << "Wrong @ Norm_Ord_Sec::cal_F1" << endl;
                    cout << "vec_3V_flag_Tz[ch_TP][index_x][index_y] : "
                         << vec_3V_flag_Tz_S.val(ch_TP, index_x, index_y)
                         << endl;
                    cout << "ch_TP : " << ch_TP << "\t index_x : " << index_x
                         << "\t index_y : " << index_y << endl;
                    cout << "\t \tJ2\t T2z\t J2_ah1\t a\t h1\t h2\t" << endl;
                    cout << "\t \t" << J2 << "\t " << T2z << "\t " << J2_ah1
                         << "\t " << a << "\t " << h1 << "\t " << h2 << endl;
                    cout << "\t \t" << J2 << "\t " << T2z << "\t " << J2_pp1
                         << "\t " << p << "\t " << p1 << "\t " << p2 << endl;
                    exit(0);
                  }
                  // double val_3V = vec_3V_Tz[ch_TP][index_x][index_y];
                  double val_3V = vec_3V_Tz_S.val(ch_TP, index_x, index_y);
                  res += (J2 + 1) * val_3V * val_3V / E_dom;
                }
              }
            }
          }
        }
      }
    }
  }
  double result = res * phase / (12.0 * (j2_a + 1.0));
  return result;
}

double Norm_Ord_Sec::cal_F2(int a) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;

  for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
    int e_p1 = SP.state_Tz[p1].E;
    int j2_p1 = SP.state_Tz[p1].j_2;
    int t2z_p1 = SP.state_Tz[p1].tz_2;
    int J2_ap1_min = std::abs(j2_a - j2_p1);
    int J2_ap1_max = j2_a + j2_p1;
    for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
      int e_p2 = SP.state_Tz[p2].E;
      int j2_p2 = SP.state_Tz[p2].j_2;
      int t2z_p2 = SP.state_Tz[p2].tz_2;
      if (pauli_exc_Tz(a, p1, p2) < 0)
        continue;
      int T2z = t2z_a + t2z_p1 + t2z_p2;

      for (int h = 0; h <= sp_core_Tz; h++) {
        int e_h = SP.state_Tz[h].E;
        int j2_h = SP.state_Tz[h].j_2;
        int t2z_h = SP.state_Tz[h].tz_2;
        for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
          int e_h1 = SP.state_Tz[h1].E;
          int j2_h1 = SP.state_Tz[h1].j_2;
          int t2z_h1 = SP.state_Tz[h1].tz_2;
          int J2_hh1_min = std::abs(j2_h - j2_h1);
          int J2_hh1_max = j2_h + j2_h1;
          for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
            int e_h2 = SP.state_Tz[h2].E;
            int j2_h2 = SP.state_Tz[h2].j_2;
            int t2z_h2 = SP.state_Tz[h2].tz_2;
            if (pauli_exc_Tz(h, h1, h2) < 0)
              continue;
            int T2z_x = t2z_h + t2z_h1 + t2z_h2;
            if (T2z != T2z_x)
              continue;
            double E_dom =
                this->hbar_omega * (e_h + e_h1 + e_h2 - e_a - e_p1 - e_p2);
            if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
              continue;
            for (int J2_ap1 = J2_ap1_min; J2_ap1 <= J2_ap1_max;
                 J2_ap1 = J2_ap1 + 2) {
              int J2_max = J2_ap1 + j2_p2;
              int J2_min = std::abs(J2_ap1 - j2_p2);
              for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
                for (int J2_hh1 = J2_hh1_min; J2_hh1 <= J2_hh1_max;
                     J2_hh1 = J2_hh1 + 2) {
                  if (Tri(J2, J2_hh1, j2_h2) < 0)
                    continue;
                  int ch_TP = TP.FindC_Tz(J2, T2z);
                  int index_x = TP.FindTP_Tz(J2, T2z, J2_ap1, a, p1, p2);
                  int index_y = TP.FindTP_Tz(J2, T2z, J2_hh1, h, h1, h2);
                  if (vec_3V_flag_Tz_S.val(ch_TP, index_y, index_x) != 1) {
                    cout << "Wrong @ Norm_Ord_Sec::cal_F2" << endl;
                    cout << "vec_3V_flag_Tz[ch_TP][index_y][index_x] : "
                         << vec_3V_flag_Tz_S.val(ch_TP, index_y, index_x)
                         << endl;
                    cout << "ch_TP : " << ch_TP << "\t index_y : " << index_y
                         << "\t index_x : " << index_x << endl;
                    cout << "\t \tJ2\t T2z\t J2_ah1\t a\t b\t c\t" << endl;
                    cout << "\t \t" << J2 << "\t " << T2z << "\t " << J2_ap1
                         << "\t " << a << "\t " << p1 << "\t " << p2 << endl;
                    cout << "\t \t" << J2 << "\t " << T2z << "\t " << J2_hh1
                         << "\t " << h << "\t " << h1 << "\t " << h2 << endl;
                    exit(0);
                  }
                  // double val_3V = vec_3V_Tz[ch_TP][index_y][index_x];
                  double val_3V = vec_3V_Tz_S.val(ch_TP, index_y, index_x);
                  res += (J2 + 1) * val_3V * val_3V / E_dom;
                }
              }
            }
          }
        }
      }
    }
  }
  double result = res * phase / (12.0 * (j2_a + 1.0));
  return result;
}

void Norm_Ord_Sec::print_U(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);
  out_file << "i\t j\t vec_U_flag\t vec_U\t vec_UT_flag\t vec_UT " << endl;
  int sp_size = SP.state_Tz.size();
  for (int i = 0; i < sp_size; i++) {
    for (int j = 0; j < sp_size; j++) {
      if (vec_U_flag[i][j] > 0 || vec_UT_flag[i][j] > 0)
        out_file << i << "\t " << j << "\t " << vec_U_flag[i][j] << "\t "
                 << vec_U[i][j] << "\t " << vec_UT_flag[i][j] << "\t"
                 << vec_UT[i][j] << endl;
    }
  }
  out_file.close();
}

void Norm_Ord_Sec::print_U_all(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);
  out_file << "index\t n\t l\t j2\t tz*2\t " << endl;
  vector<State_SP> state_Tz;
  state_Tz = SP.state_Tz;
  std::stable_sort(
      state_Tz.begin(), state_Tz.end(), [](State_SP a, State_SP b) {
        return (100000 * a.E + 10000 * a.n - 1000 * a.l - 100 * a.j_2 +
                a.tz_2) <
               (100000 * b.E + 10000 * b.n - 1000 * b.l - 100 * b.j_2 + b.tz_2);
      });
  int size_sp = state_Tz.size();
  vector<int> this2other(size_sp);
  for (int i = 0; i < size_sp; i++) {
    int flag = -1;
    int n = SP.state_Tz[i].n;
    int l = SP.state_Tz[i].l;
    int j2 = SP.state_Tz[i].j_2;
    int tz2 = SP.state_Tz[i].tz_2;
    for (int j = 0; j < size_sp; j++) {
      int n_t = state_Tz[j].n;
      int l_t = state_Tz[j].l;
      int j2_t = state_Tz[j].j_2;
      int tz2_t = state_Tz[j].tz_2;
      if (n_t != n || l_t != l)
        continue;
      if (j2_t != j2 || tz2_t != tz2)
        continue;
      flag = j;
      break;
    }
    if (flag < 0) {
      cout << "Wrong happened @Norm_Ord_Sec::print_2V_all" << endl;
      exit(0);
    }
    this2other[i] = flag + 1;
  }

  for (int i = 0; i < size_sp; i++) {
    int index = i + 1;
    int n = state_Tz[i].n;
    int l = state_Tz[i].l;
    int j2 = state_Tz[i].j_2;
    int tz_2 = state_Tz[i].tz_2;
    out_file << index << "\t " << n << "\t " << l << "\t " << j2 << "\t "
             << tz_2 << endl;
  }
  out_file << "--++++++++++++++++++++++++++++++++++++++++++++++++++--" << endl;
  out_file << "i\t j\t vec_U_flag\t vec_U\t vec_UT_flag\t vec_UT " << endl;
  int sp_size = SP.state_Tz.size();
  for (int i = 0; i < sp_size; i++) {
    for (int j = 0; j < sp_size; j++) {
      if (vec_U_flag[i][j] > 0 || vec_UT_flag[i][j] > 0)
        out_file << this2other[i] << "\t " << this2other[j] << "\t "
                 << vec_U_flag[i][j] << "\t " << vec_U[i][j] << "\t "
                 << vec_UT_flag[i][j] << "\t" << vec_UT[i][j] << endl;
    }
  }
  out_file.close();
}

void Norm_Ord_Sec::print_2V(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);

  int ch_size = DP.state_Tz.size();

  for (int ch = 0; ch < ch_size; ch++) {
    out_file << " --- DP channel_Tz : " << ch << " --- " << endl;

    out_file << "i\t j\t vec_2V2N_flag\t vec_2V2N\t vec_2V3N_flag\t vec_2V3N "
             << endl;
    int size = DP.state_Tz[ch].size();
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (vec_2V3N_flag[ch][i][j] > 0 || vec_2V2N_flag[ch][i][j] > 0)
          out_file << i << "\t " << j << "\t " << vec_2V2N_flag[ch][i][j]
                   << "\t " << vec_2V2N[ch][i][j] << "\t "
                   << vec_2V3N_flag[ch][i][j] << "\t" << vec_2V3N[ch][i][j]
                   << endl;
      }
    }
  }

  out_file.close();
}

void Norm_Ord_Sec::print_2V_all(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);
  out_file << "index\t n\t l\t j2\t tz*2\t " << endl;

  vector<State_SP> state_Tz;
  state_Tz = SP.state_Tz;
  std::stable_sort(
      state_Tz.begin(), state_Tz.end(), [](State_SP a, State_SP b) {
        return (100000 * a.E + 10000 * a.n - 1000 * a.l - 100 * a.j_2 +
                a.tz_2) <
               (100000 * b.E + 10000 * b.n - 1000 * b.l - 100 * b.j_2 + b.tz_2);
      });
  int size_sp = state_Tz.size();
  vector<int> this2other(size_sp);
  for (int i = 0; i < size_sp; i++) {

    int n = SP.state_Tz[i].n;
    int l = SP.state_Tz[i].l;
    int j2 = SP.state_Tz[i].j_2;
    int tz2 = SP.state_Tz[i].tz_2;
    int flag = -1;
    for (int j = 0; j < size_sp; j++) {
      int n_t = state_Tz[j].n;
      int l_t = state_Tz[j].l;
      int j2_t = state_Tz[j].j_2;
      int tz2_t = state_Tz[j].tz_2;
      if (n_t != n || l_t != l)
        continue;
      if (j2_t != j2 || tz2_t != tz2)
        continue;
      flag = j;
      break;
    }
    if (flag < 0) {
      cout << "Wrong happened @Norm_Ord_Sec::print_2V_all" << endl;
      exit(0);
    }
    this2other[i] = flag + 1;
    cout << "this2other[" << i << "] = " << this2other[i] << endl;
  }

  for (int i = 0; i < size_sp; i++) {
    int index = i + 1;
    int n = state_Tz[i].n;
    int l = state_Tz[i].l;
    int j2 = state_Tz[i].j_2;
    int tz_2 = state_Tz[i].tz_2;
    out_file << index << "\t " << n << "\t " << l << "\t " << j2 << "\t "
             << tz_2 << endl;
  }
  out_file << "--++++++++++++++++++++++++++++++++++++++++++++++++++--" << endl;

  out_file << " T2z\t J2\t Par\t\t a\t b\t c\t d\t\t vec_2V2N\t vec_2V2N_unN\t"
              " vec_2V3N "
           << endl;
  int ch_size = DP.state_Tz.size();
  for (int ch = 0; ch < ch_size; ch++) {
    // out_file << " --- DP channel_Tz : " << ch << " --- " << endl;
    int size = DP.state_Tz[ch].size();
    for (int i = 0; i < size; i++) {
      int a = DP.state_Tz[ch][i].a;
      int b = DP.state_Tz[ch][i].b;
      int a_t = this2other[a];
      int b_t = this2other[b];
      int J2 = DP.state_Tz[ch][i].J_2;
      int T2z = DP.state_Tz[ch][i].T_2z;
      int Par = DP.state_Tz[ch][i].Par;
      for (int j = 0; j < size; j++) {
        int c = DP.state_Tz[ch][j].a;
        int d = DP.state_Tz[ch][j].b;
        int c_t = this2other[c];
        int d_t = this2other[d];
        if (vec_2V3N_flag[ch][i][j] > 0 || vec_2V2N_flag[ch][i][j] > 0)
          out_file << T2z << "\t " << J2 << "\t " << Par << "\t\t " << a_t
                   << "\t " << b_t << "\t " << c_t << "\t " << d_t << "\t   "
                   << std::setw(12) << vec_2V2N[ch][i][j] << "\t \t"
                   << std::setw(12) << vec_2V2N_unnorm[ch][i][j] << "\t \t"
                   << std::setw(12) << vec_2V3N[ch][i][j] << endl;
      }
    }
  }

  out_file.close();
}

void Norm_Ord_Sec::print_2V_Luigi(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);
  out_file << "n\t l\t j2\t index " << endl;

  vector<State_SP> state;
  state = SP.state;
  std::stable_sort(state.begin(), state.end(), [](State_SP a, State_SP b) {
    return (100000 * a.E + 10000 * a.n - 1000 * a.l - 100 * a.j_2) <
           (100000 * b.E + 10000 * b.n - 1000 * b.l - 100 * b.j_2);
  });
  int size_sp = SP.state_Tz.size();
  int size_sp_t = state.size();
  vector<int> this2other(size_sp);
  for (int i = 0; i < size_sp; i++) {

    int n = SP.state_Tz[i].n;
    int l = SP.state_Tz[i].l;
    int j2 = SP.state_Tz[i].j_2;
    // int tz2 = SP.state_Tz[i].tz_2;
    int flag = -1;
    for (int j = 0; j < size_sp_t; j++) {
      int n_t = state[j].n;
      int l_t = state[j].l;
      int j2_t = state[j].j_2;
      //  int tz2_t = state_Tz[j].tz_2;
      if (n_t != n || l_t != l)
        continue;
      if (j2_t != j2)
        continue;
      flag = j;
      break;
    }
    if (flag < 0) {
      cout << "Wrong happened @Norm_Ord_Sec::print_2V_all" << endl;
      exit(0);
    }
    this2other[i] = flag + 1;
    cout << "this2other[" << i << "] = " << this2other[i] << endl;
  }

  for (int i = 0; i < size_sp_t; i++) {
    int index = i + 1;
    int n = state[i].n;
    int l = state[i].l;
    int j2 = state[i].j_2;
    // int tz_2 = state_Tz[i].tz_2;
    out_file << n << "\t " << l << "\t " << j2 << "\t " << index << endl;
  }
  out_file << "--++++++++++++++++++++++++++++++++++++++++++++++++++--" << endl;
  out_file << "\t <ab|V|cd> Symmetry : a<=b; c<=d; a<=c" << endl;
  out_file << "Type 1 : <nn|V|nn>;  Type 2 : <pn|V|pn>;  Type 3 : <pn|V|np>"
           << endl;
  out_file << endl;
  out_file
      << "Type\t T2z\t J2\t Par\t\t a\t b\t c\t d\t\t vec_2V2N\t vec_2V2N_unN\t"
         "\t vec_2V3N "
      << endl;
  int ch_size = DP.state_Tz.size();
  for (int ty = 0; ty < 3; ty++) {
    for (int Par = 0; Par <= 1; Par++) {
      for (int ch = 0; ch < ch_size; ch++) {
        // out_file << " --- DP channel_Tz : " << ch << " --- " << endl;
        int T2z = DP.channel_Tz[ch].T_2z;
        if (ty == 0 && T2z != -2)
          continue;
        if (ty == 1 && T2z != 0)
          continue;
        if (ty == 2 && T2z != 0)
          continue;
        int size = DP.state_Tz[ch].size();
        for (int i = 0; i < size; i++) {
          int a = DP.state_Tz[ch][i].a;
          int b = DP.state_Tz[ch][i].b;
          int t2z_a = SP.state_Tz[a].tz_2;
          // int t2z_b = SP.state_Tz[b].tz_2;
          int a_t = this2other[a];
          int b_t = this2other[b];
          if (b_t < a_t)
            continue;
          if (ty == 1 && t2z_a < 0)
            continue;
          if (ty == 2 && t2z_a < 0)
            continue;
          int J2 = DP.state_Tz[ch][i].J_2;
          // int T2z = DP.state_Tz[ch][i].T_2z;
          int Par_i = DP.state_Tz[ch][i].Par;
          if (Par != Par_i)
            continue;

          for (int j = 0; j < size; j++) {
            int Par_j = DP.state_Tz[ch][i].Par;
            if (Par != Par_j)
              continue;
            int c = DP.state_Tz[ch][j].a;
            int d = DP.state_Tz[ch][j].b;
            int c_t = this2other[c];
            int d_t = this2other[d];
            int t2z_c = SP.state_Tz[c].tz_2;
            // int t2z_d = SP.state_Tz[b].tz_2;
            if (ty == 1 && t2z_c < 0)
              continue;
            if (ty == 2 && t2z_c > 0)
              continue;
            if (d_t < c_t)
              continue;
            if (vec_2V3N_flag[ch][i][j] > 0 || vec_2V2N_flag[ch][i][j] > 0)
              out_file << ty + 1 << "\t " << T2z << "\t " << J2 << "\t "
                       << Par << "\t\t " << a_t << "\t "
                       << b_t << "\t " << c_t << "\t " << d_t << "\t   "
                       << std::setw(12) << vec_2V2N[ch][i][j] << "\t \t"
                       << std::setw(12) << vec_2V2N_unnorm[ch][i][j] << "\t \t"
                       << std::setw(12) << vec_2V3N[ch][i][j] << endl;
          }
        }
      }
    }
  }

  out_file.close();
}

int Norm_Ord_Sec::pauli_exc_Tz(int a, int b, int c) {
  int res = 1;
  if (a == b && b == c && SP.state_Tz[a].l == 0)
    res = -1;
  return res;
}
//==========--------------------------------------=============//

void Norm_Ord_Core::init(double hw, int sp_max, int sp_core, int sp_P_beg,
                         int sp_P_end, int sp_Q_beg, int sp_Q_end) {
  this->hbar_omega = hw;
  cout << "!!! Norm_Ord_Core Hbar_Omega = " << hw << endl;
  this->E_cut = 0; // E_cut;
  this->sp_max = sp_max;
  this->sp_core = sp_core;
  this->sp_P_beg = sp_P_beg;
  this->sp_P_end = sp_P_end;
  this->sp_Q_beg = sp_Q_beg;
  this->sp_Q_end = sp_Q_end;

  this->sp_max_Tz = sp_max * 2;
  this->sp_core_Tz = sp_core * 2 + 1;
  this->sp_P_beg_Tz = sp_P_beg * 2;
  this->sp_P_end_Tz = sp_P_end * 2 + 1;
  this->sp_Q_beg_Tz = sp_Q_beg * 2;
  this->sp_Q_end_Tz = sp_Q_end * 2 + 1;
  int size_ch = DP.state_Tz.size();
  // cout << "DP channel : " << size_ch << endl;

  vec_2V3N.resize(size_ch);
  vec_2V3N_flag.resize(size_ch);
  for (int ch = 0; ch < size_ch; ch++) {
    int size = DP.state_Tz[ch].size();
    vec_2V3N[ch].resize(size);
    vec_2V3N_flag[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2V3N[ch][i].resize(size);
      vec_2V3N_flag[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec_2V3N[ch][i][j] = 0.0;
        vec_2V3N_flag[ch][i][j] = -1;
      }
    }
  }
}
double Norm_Ord_Core::cal2f3_Tz(State_DP &T2P_i, State_DP &T2P_f) {
  double res = 0.0;

  int T2ab_z_i = T2P_i.T_2z;
  int T2ab_z_f = T2P_f.T_2z;
  int J2ab_i = T2P_i.J_2;
  int J2ab_f = T2P_f.J_2;
  if (T2ab_z_i != T2ab_z_f)
    return res;
  if (J2ab_i != J2ab_f)
    return res;
  int T2ab_z = T2ab_z_i;
  int J2ab = J2ab_i;
  int a = T2P_i.a;
  int b = T2P_i.b;
  int d = T2P_f.a;
  int e = T2P_f.b;

  int h_core = this->sp_core_Tz;

  // cout << "--- -- a = " << a << "\t b = " << b << "\t\t d = " << d
  //     << "\t e = " << e << endl;

  // double result = 0.0;
  for (int h = 0; h <= h_core; h++) {
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2max = j2_h + J2ab;
    int J2min = std::abs(j2_h - J2ab);
    int T2z = T2ab_z + t2z_h;
    for (int J2 = J2min; J2 <= J2max; J2 = J2 + 2) {
      int ch_Tz = TP.FindC_Tz(J2, T2z);
      /*cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " <<
      J2ab
           << endl;
      cout << "a = " << a << "\t b = " << b << "\t h = " << h <<
      endl;/**/
      if (SP.state_Tz[a].l == 0 && a == b && b == h)
        continue;
      if (SP.state_Tz[d].l == 0 && d == e && e == h)
        continue; /**/
      int x = TP.FindTP_Tz(J2, T2z, J2ab, a, b, h);
      // cout << "d = " << d << "\t e = " << e << "\t h = " << h <<
      // endl;/**/
      int y = TP.FindTP_Tz(J2, T2z, J2ab, d, e, h);
      if (x < 0) {
        cout << "Wrong at cal2f3_Tz x< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "a = " << a << "\t b = " << b << "\t h = " << h << endl;
        exit(0);
      }
      if (y < 0) {
        cout << "Wrong at cal2f3_Tz y< 0" << endl;
        cout << "J2 = " << J2 << "\t T2z = " << T2z << "\t J2ab = " << J2ab
             << endl;
        cout << "d = " << d << "\t e = " << e << "\t h = " << h << endl;
        exit(0);
      }
      // cout << "ch_Tz = " << ch_Tz << "\t x = " << x << "\t y = " << y
      // << endl;
      // if (this->vec_3V_flag_Tz[ch_Tz][x][y] != 1) {
      if (this->vec_3V_flag_Tz_S.val(ch_Tz, x, y) != 1) {
        cout << "Wrong happened @ Norm_Ord::cal_val2f3_Tz 1" << endl;
        cout << "this->vec_flag_Tz[ch_Tz][x][y] != 1" << endl;
        exit(0);
      }
      // double mat_t = this->vec_3V_Tz[ch_Tz][x][y];
      double mat_t = this->vec_3V_Tz_S.val(ch_Tz, x, y);
      // cout << "mat_t = " << mat_t << endl;
      double mat = (J2 + 1.0) / (J2ab + 1.0) * mat_t;
      res += mat;
    }
  }
  return res;
}
void Norm_Ord_Core::build_P() {
  // 1: <ab(h0)|V3N|cd(h0)>

  for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
    int t2z_a = SP.state_Tz[a].tz_2;
    int j2_a = SP.state_Tz[a].j_2;
    int l_a = SP.state_Tz[a].l;
    int e_a = SP.state_Tz[a].E;

    for (int b = this->sp_P_beg_Tz; b <= this->sp_P_end_Tz; b++) {
      int t2z_b = SP.state_Tz[b].tz_2;
      int j2_b = SP.state_Tz[b].j_2;
      int l_b = SP.state_Tz[b].l;
      int e_b = SP.state_Tz[b].E;

      int T2z_x = t2z_a + t2z_b;
      int Par_x = (l_a + l_b) % 2;
      int J2_min = std::abs(j2_a - j2_b);
      int J2_max = j2_a + j2_b;

      for (int c = this->sp_P_beg_Tz; c <= this->sp_P_end_Tz; c++) {
        int t2z_c = SP.state_Tz[c].tz_2;
        int j2_c = SP.state_Tz[c].j_2;
        int l_c = SP.state_Tz[c].l;
        int e_c = SP.state_Tz[c].E;
        for (int d = this->sp_P_beg_Tz; d <= this->sp_P_end_Tz; d++) {
          int t2z_d = SP.state_Tz[d].tz_2;
          int j2_d = SP.state_Tz[d].j_2;
          int l_d = SP.state_Tz[d].l;
          int e_d = SP.state_Tz[d].E;
          // cout << "\t a : " << a << "\t b : " << b << "\t c : " << c
          //      << "\t d : " << endl;
          if (T2z_x != (t2z_c + t2z_d))
            continue;
          if (Par_x != (l_c + l_d) % 2)
            continue;
          // int E_dom_int = std::abs(e_a + e_h - e_p1 - e_p2);
          // if (E_dom_int > this->E_cut)
          //   continue;
          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_c, j2_d) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, a, b);
            int index_y = DP.FindDP_Tz(J2, T2z_x, c, d);

            double mat_3 = this->cal2f3_Tz(DP.state_Tz[ch_Tz][index_x],
                                           DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;
          }
        }
      }
    }
  } /**/
}
void Norm_Ord_Core::build_Core() {
  // Totally U has two piece
  // 1: <h1a(h0)|V3N|bc(h0)>  2: <h1h2(h0)|V3N|bc(h0)>
  // int size_ch = DP.state_Tz.size();
  // // cout << "DP channel : " << size_ch << endl;
  //
  // vec_2V3N.resize(size_ch);
  // vec_2V3N_flag.resize(size_ch);
  // for (int ch = 0; ch < size_ch; ch++) {
  //   int size = DP.state_Tz[ch].size();
  //   vec_2V3N[ch].resize(size);
  //   vec_2V3N_flag[ch].resize(size);
  //   for (int i = 0; i < size; i++) {
  //     vec_2V3N[ch][i].resize(size);
  //     vec_2V3N_flag[ch][i].resize(size);
  //     for (int j = 0; j < size; j++) {
  //       vec_2V3N[ch][i][j] = 0.0;
  //       vec_2V3N_flag[ch][i][j] = -1;
  //     }
  //   }
  // }
  // 1: <h1a(h0)|V3N|bc(h0)> & <ah1(h0)|V3N|bc(h0)>

  for (int h1 = 0; h1 <= this->sp_core_Tz; h1++) {
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    int j2_h1 = SP.state_Tz[h1].j_2;
    int l_h1 = SP.state_Tz[h1].l;
    int e_h1 = SP.state_Tz[h1].E;
    for (int a = this->sp_P_beg_Tz; a <= this->sp_P_end_Tz; a++) {
      int t2z_a = SP.state_Tz[a].tz_2;
      int j2_a = SP.state_Tz[a].j_2;
      int l_a = SP.state_Tz[a].l;
      int e_a = SP.state_Tz[a].E;

      int T2z_x = t2z_a + t2z_h1;
      int Par_x = (l_a + l_h1) % 2;
      int J2_min = std::abs(j2_a - j2_h1);
      int J2_max = j2_a + j2_h1;

      for (int c = this->sp_P_beg_Tz; c <= this->sp_P_end_Tz; c++) {
        int t2z_c = SP.state_Tz[c].tz_2;
        int j2_c = SP.state_Tz[c].j_2;
        int l_c = SP.state_Tz[c].l;
        int e_c = SP.state_Tz[c].E;
        for (int d = this->sp_P_beg_Tz; d <= this->sp_P_end_Tz; d++) {
          int t2z_d = SP.state_Tz[d].tz_2;
          int j2_d = SP.state_Tz[d].j_2;
          int l_d = SP.state_Tz[d].l;
          int e_d = SP.state_Tz[d].E;
          if (T2z_x != (t2z_c + t2z_d))
            continue;
          if (Par_x != (l_c + l_d) % 2)
            continue;
          // int E_dom_int = std::abs(e_a + e_h - e_p1 - e_p2);
          // if (E_dom_int > this->E_cut)
          //   continue;
          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_c, j2_d) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, h1, a);
            int index_y = DP.FindDP_Tz(J2, T2z_x, c, d);

            double mat_3 = this->cal2f3_Tz(DP.state_Tz[ch_Tz][index_x],
                                           DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;

            index_x = DP.FindDP_Tz(J2, T2z_x, a, h1);
            mat_3 = this->cal2f3_Tz(DP.state_Tz[ch_Tz][index_x],
                                    DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;
          }
        }
      }
    }
  } /**/
  // 2: <h1h2(h0)|V3N|bc(h0)>
  for (int h1 = 0; h1 <= this->sp_core_Tz; h1++) {
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    int j2_h1 = SP.state_Tz[h1].j_2;
    int l_h1 = SP.state_Tz[h1].l;
    int e_h1 = SP.state_Tz[h1].E;
    for (int h2 = 0; h2 <= this->sp_core_Tz; h2++) {
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      int j2_h2 = SP.state_Tz[h2].j_2;
      int l_h2 = SP.state_Tz[h2].l;
      int e_h2 = SP.state_Tz[h2].E;

      int T2z_x = t2z_h1 + t2z_h2;
      int Par_x = (l_h1 + l_h2) % 2;
      int J2_min = std::abs(j2_h1 - j2_h2);
      int J2_max = j2_h1 + j2_h2;

      for (int c = this->sp_P_beg_Tz; c <= this->sp_P_end_Tz; c++) {
        int t2z_c = SP.state_Tz[c].tz_2;
        int j2_c = SP.state_Tz[c].j_2;
        int l_c = SP.state_Tz[c].l;
        int e_c = SP.state_Tz[c].E;
        for (int d = this->sp_P_beg_Tz; d <= this->sp_P_end_Tz; d++) {
          int t2z_d = SP.state_Tz[d].tz_2;
          int j2_d = SP.state_Tz[d].j_2;
          int l_d = SP.state_Tz[d].l;
          int e_d = SP.state_Tz[d].E;

          if (T2z_x != (t2z_c + t2z_d))
            continue;
          // int E_dom_int = std::abs(e_h1 + e_h2 - e_a - e_p);
          // if (E_dom_int > this->E_cut)
          //  continue;
          // if (Par_x != (l_a + l_p) % 2)
          //   continue;

          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_c, j2_d) < 0)
              continue;

            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            int index_x = DP.FindDP_Tz(J2, T2z_x, h1, h2);

            int index_y = DP.FindDP_Tz(J2, T2z_x, c, d);

            double mat_3 = this->cal2f3_Tz(DP.state_Tz[ch_Tz][index_x],
                                           DP.state_Tz[ch_Tz][index_y]);
            vec_2V3N[ch_Tz][index_x][index_y] = mat_3;
            vec_2V3N_flag[ch_Tz][index_x][index_y] = 1;
            vec_2V3N[ch_Tz][index_y][index_x] = mat_3;
            vec_2V3N_flag[ch_Tz][index_y][index_x] = 1;
          }
        }
      }
    }
  } /**/
}
void Norm_Ord_Core::print_norm_core(std::string filename, int h_core_Tz,
                                    int p_beg_Tz, int p_end_Tz, double hw,
                                    int reg_pow, int Lambda_Cut,
                                    vector<int> oc_o2z) {
  // y= oc_o2z[x] : y is state_Tz;  x is from others
  vector<int> oc_z2o;
  int size_oc = oc_o2z.size();
  oc_z2o.resize(size_oc);
  for (int x = 0; x < size_oc; x++) {
    int y = oc_o2z[x];
    oc_z2o[y] = x;
  }

  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t gA = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hw << endl;
  file << "\t The Power of Regulator : " << reg_pow << endl;
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t";
  for (int i = p_end; i >= p_beg; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    int y = oc_o2z[i];
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[y].n << "\t";
    file << SP.state[y].l << "\t";
    file << SP.state[y].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "\t Num: " << p_end_Tz - 0 + 1 << endl;
  file << "\tindex\t n\t l\t 2*j\t 2*tz\t E (MeV)" << endl;
  for (int i = p_end_Tz; i >= 0; i--) {
    int y = std::floor(i / 2.0);
    int x = oc_z2o[y];
    file << "\t" << x + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    // file << std::setw(6) << vec_1f3_C[i] << endl;
    file << std::setw(6) << 0 << endl;
  }
  file << endl;
  file << " === Normal-Ordering Two-Body ME : sum over {core} orbits=== "
       << endl;
  file << "\tMT=-1 : <nn|V|nn>; MT=1 <pp|V|pp>; MT=0: <pn|V|pn>" << endl;
  file << "\t<(a'b')Jab mta' mtb' MT |V| (ab)Jab mta mtb MT>" << endl;
  file << endl;
  file << "a'\tb'\ta\tb\tJab\t MT\t ME (MeV)" << endl;
  SP.print_Tz();
  for (int a = 0; a <= this->sp_P_end_Tz; a++) {
    int t2z_a = SP.state_Tz[a].tz_2;
    int j2_a = SP.state_Tz[a].j_2;
    int l_a = SP.state_Tz[a].l;
    int e_a = SP.state_Tz[a].E;
    // cout << "\t \t a : " << a << "\t tz = " << t2z_a << endl;
    for (int b = 0; b <= this->sp_P_end_Tz; b++) {
      int t2z_b = SP.state_Tz[b].tz_2;
      int j2_b = SP.state_Tz[b].j_2;
      int l_b = SP.state_Tz[b].l;
      int e_b = SP.state_Tz[b].E;

      int T2z_x = t2z_a + t2z_b;
      int Par_x = (l_a + l_b) % 2;
      int J2_min = std::abs(j2_a - j2_b);
      int J2_max = j2_a + j2_b;

      for (int c = 0; c <= this->sp_P_end_Tz; c++) {
        int t2z_c = SP.state_Tz[c].tz_2;
        int j2_c = SP.state_Tz[c].j_2;
        int l_c = SP.state_Tz[c].l;
        int e_c = SP.state_Tz[c].E;
        for (int d = 0; d <= this->sp_P_end_Tz; d++) {
          int t2z_d = SP.state_Tz[d].tz_2;
          int j2_d = SP.state_Tz[d].j_2;
          int l_d = SP.state_Tz[d].l;
          int e_d = SP.state_Tz[d].E;

          if (T2z_x != (t2z_c + t2z_d))
            continue;
          if (Par_x != (l_c + l_d) % 2)
            continue;
          // int E_dom_int = std::abs(e_a + e_h - e_p1 - e_p2);
          // if (E_dom_int > this->E_cut)
          //   continue;
          for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
            if (this->Tri(J2, j2_c, j2_d) < 0)
              continue;
            int ch_Tz = DP.FindC_Tz(J2, T2z_x);
            /*cout << "J2 : " << J2 << "\t T2z_x : " << T2z_x
                 << "\t ch_Tz : " << ch_Tz << "\t a : " << a << "\t b : " << b
                 << "\t c : " << c << "\t d : " << d << endl;/**/
            /*cout << "\t t2z : a : " << t2z_a << "\t b : " << t2z_b
                 << "\t c : " << t2z_c << "\t d : " << t2z_d << endl;/**/
            int index_x = DP.FindDP_Tz(J2, T2z_x, a, b);

            int index_y = DP.FindDP_Tz(J2, T2z_x, c, d);

            int ch = ch_Tz;
            int i = index_x;
            int j = index_y;
            if (vec_2V3N_flag[ch][i][j] < 0)
              continue;
            int a_x = floor(a / 2);
            int b_x = floor(b / 2);
            int c_x = floor(c / 2);
            int d_x = floor(d / 2);
            /* Only considering <pn|V|pn> */
            if (T2z_x == 0 && t2z_a < 0)
              continue;
            if (T2z_x == 0 && t2z_c < 0)
              continue;
            /*file << ap + 1 << "\t" << bp + 1 << "\t" << a + 1 << "\t" << b +
               1
                 << "\t" << J2ab / 2 << "\t" << T2ab_z / 2 << "\t";/**/
            file << oc_z2o[a_x] + 1 << "\t" << oc_z2o[b_x] + 1 << "\t"
                 << oc_z2o[c_x] + 1 << "\t" << oc_z2o[d_x] + 1 << "\t" << J2 / 2
                 << "\t" << T2z_x / 2 << "\t";
            file << vec_2V3N[ch][i][j];
            // file << "\t \t " << vec_2V3N_flag[ch][i][j];
            file << endl;
          }
        }
      }
    }
  } /**/

  file.close();
} /**/
int Norm_Ord_Core::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}
